数值分析

基础知识

实数的浮点表示

机器数的格式

浮点数的表示：详见二进制转换成其它格式数据，此处的浮点数精度为64位，即指数部分为11位，尾数部分为52位。标准化的IEEE浮点数表示为：

\pm 1. b b \dots b \times 2^{p}

$\varepsilon_{mach}=2^{-52}$

IEEE舍入最近法则

如果在二进制数右边的第53位是0，则向下舍去（在第52位后面截断）。
如果第53位是1，则向上进位（在第52位上加1)；但是如果在53位的右边的所有已知位都是0，此时当且仅当第52位是1时才在第52位上加1，否则在52位后面截断。

定义：将IEEE双精度浮点数字记做x，利用舍入最近法则记做fl(x)。

举例 $\left( {1001.\overline {0110} } \right)_2$ ，计算舍入误差。

$1.001\overline {0110} \times {2^3} = 1.0010\overline {1100} \times {2^3}$ $1 \times {2^{ - 52}} \times {2^3}$ $0.\overline {1100} \times {2^{ - 52}} \times {2^3} = 0.\overline {1100} \times {2^{ - 49}}$ ，

$0.\overline {1100}$ ，可以这样计算：

\begin{aligned} x = 0. \overset{―}{1100} \\ 2^{4} x = 1100. \overset{―}{1100} \\ (2^{4} - 1) x = 1100 \Rightarrow x = 0.8 \end{aligned}

所以舍入误差为

f l (x) - x = - 0.8 \times 2^{- 49} + 1 \times 2^{- 52} \times 2^{3} = 0.2 \times 2^{- 49}

相对舍入误差

在IEEE机器算术模型中，fl(x)的相对舍入误差不会比机器精度的一半大：
$\frac{| f l (x) - x |}{| x |} \leq \frac{1}{2} ε_{m a c h}$

机器数的表示

举例：fl(9.4)对应的存储在机器中的表达？

${\left( {10000000010} \right)_2}$ ，再加上符号位0，则为

{(010000000010)}_{2} = {(402)}_{16}

${\left( {{{2CCCCCCCCCCCD}}} \right)_{16}}$ 。

整体16进制表示为4022CCCCCCCCCCCD。

特殊数的十六进制表示

机器数	例子	十六进制数
+Inf	1/0	7FF0000000000000
-Inf	-1/0	FFF0000000000000
NaN	0/0	FFFxxxxxxxxxxxxx
+0		0000000000000000
-0		8000000000000000

其中x表示不全为0的位，+0和-0在计算中相同，只有符号位不同。

${2^{ - 52}} \times {2^{ - 1022}} = {2^{ - 1074}}$ ，这是最小可表达的数字，在此之下的可表达双精度数字完全不能被表示。

浮点数加法

$2^{-53}$ $2^{-53}$ ，虽然53位是1，但是右边全是0且第52位是0，浮点数表示为0，即

1 + 2^{- 53} = 1

举例 $1+3\times2^{-53}$ )-1

由于求和后的53位是1，又52位也是1，所以需要在第52位上加上1，即最后结果为：

$2^{-51}$ ，因为舍入误差导致最后的结果不为0。

有效数字缺失

举例 $x^2+9^{12}x=3$ 的两个根。

$x = {{ - b \pm \sqrt {{b^2} - 4ac} } \over {2a}}$ ，负号的值可以找到，但是对于正号的值，存在舍入误差会为0，这就是有效数字缺失（两个近似的数字相减）。可以将正号时的求根公式转成

\begin{aligned} x_{2} = \frac{- b + \sqrt{b^{2} - 4 a c}}{2 a} = \frac{(- b + \sqrt{b^{2} - 4 a c}) (- b - \sqrt{b^{2} - 4 a c})}{2 a (- b - \sqrt{b^{2} - 4 a c})} \\ = \frac{2 c}{- b - \sqrt{b^{2} - 4 a c}} \end{aligned}

$x_2=1.062\times10^{-11}$ 。

求解方程

二分法

二分法

给定初始区间[a,b] 使得f(a)f(b)<0 while (b-a)/2 > TOL c = (a+b) / 2 if f(c) = 0, stop,end if f(a)f(c) < 0 b=c else a=c end end 最终的区间[a,b]中包含一个根近似根为(a+b)/2

二分法的求解误差

| x_{c} - r | < \frac{b - a}{2^{n + 1}}

$x_c=(a_n+b_n)/2$ 是解的最优估计值。

定义 $0.5\times 10^{-p}$ ，解精确到小数点后p位。

举例：使用二分法求解函数在区间[0,1]上的解，要求解精确到小数点后6位，则至少要迭代多少次。

\frac{b - a}{2^{n + 1}} = \frac{1}{2^{n + 1}} < 0.5 \times 10^{- 6} \Rightarrow n \geq 20

至少要20次迭代。

不动点迭代

函数的不动点

定义：当g(r)=r，实数r是函数g的不动点。

不动点迭代

$x_0$ =初始估计
$x_{i+1}=g(x_i)$ ，其中i=0,1,2,...

举例 $x^3+x-1=0$ 的根。

$x=1-x^3$ $g(x)=1-x^3$ $g(x) = \root 3 \of {1 - x} = x$ 。

$2x^3$ ，可以得到

\begin{aligned} 3 x^{3} + x = 1 + 2 x^{3} \\ x (1 + 3 x^{2}) = 1 + 2 x^{3} \\ x = \frac{1 + 2 x^{3}}{1 + 3 x^{2}} = g (x) \end{aligned}

$x_0=0.5$ 作为初始值，但是第一种没有成功收敛（在0和1之间来回跳跃），第二种和第三种成功收敛并且第三种收敛得更快。

不动点迭代几何

使用cobweb图描述不动点迭代的几何图示

可以看到对于g(x)=1-x^3而言，在迭代时会出现发散的情况，因此最后无法收敛，而是否能够收敛与函数在不动点附近的斜率有关。

不动点迭代的线性收敛

定义 $e_i$ 表示迭代过程中第i步时的误差，如果

lim_{i \to \infty} \frac{e_{i + 1}}{e_{i}} = S < 1

该方法被称为满足线性收敛，收敛速度为S。

定理 $|g'(r)|$ <1，则不动点迭代对于一个足够接近r的初始估计，以速度S线性收敛到不动点r。

$x_i$ $x_i$ $c_i$ ，满足

g (x_{i}) - g (r) = g^{'} (c_{i}) (x_{i} - r)

$x_{i+1}=g(x_i)$ $r=g(r)$ 可得

x_{i + 1} - r = g^{'} (c_{i}) (x_{i} - r)

$e_i=|x_i-r|$ ，则有

e_{i + 1} = | g^{'} (c_{i}) | e_{i}

$x_i$ $c_i$ 也在该区间，因而

e_{i + 1} \leq \frac{S + 1}{2} e_{i}

$\mathop {\lim }\limits_{i \to \infty } {x_i} = r$ 。

定义：如果迭代方法对于一个足够接近r的初值能收敛到r，该迭代方法被称为局部收敛到r。

举例：解释为什么不动点迭代g(x)=cos(x)收敛

g'(r)=-sin(r)=-sin(0.74)≈-0.67，绝对值比1小，对应解r≈0.74会把附近的初始猜测值吸引过来。

举例：使用不动点迭代找出方程 cos x = sin x 的根。

在方程两边加上x可得

x + \cos x - \sin x = x

定义g(x)=x + cos x - sin x，可用不动点迭代法求得x=0.7853981633957986。

精度的极限

前向和后向误差

举例 $f(x)=x^3-2x^2+{4\over3}x-{8\over 27}$ 的根，精确到小数点后6位。

在使用二分法时，在16步计算后便出现错误，f(0.6666641)=0，无法精确到小数点后6位。即使使用MATLAB提供的fzero.m函数计算仍然无法找到5位以上的正确数位。

定义 $x_a$ $x_a$ 后向误差 $|f(x_a)|$ 前向误差 $|r-x_a|$ 。

$x_a$ $|f(x_a)|$ $|r-x_a|$ 。

${\varepsilon _{mach}} \approx 2.2 \times {10^{ - 16}}$ $10^{-5}$ ，双精度数不能在机器精度的相对误差以下可靠计算，因为后向误差不能被可靠降低，同时前向误差也不能减小。

方程求解方法的终止条件可以基于前向或者后向误差。还有其他相关的终止条件，诸如计算时间的上限问题的上下文对于我们的选择具有指导作用。

威尔金森多项式

由于函数在多重根位置上f'为0，所以在重根附近形状很平。正因为如此，分离重根可能会遇到困难，但是多重根问题仅仅是冰山一角，没有多重根问题时也可能也会出现问题，比如威尔金森多项式。

W (x) = (x - 1) (x - 2) \dots (x - 20)

根搜索的敏感性

如果在输人中是一个小误差，f在这种情况下对问题进行求解，造成输出中的大问题，这样的问题被称为敏感性问题。

假设对输入做了一个小变化εg(x)，其中ε很小，令Δr对应根中的变化，则

f (r + Δ r) + ε g (r + Δ r) = 0

进行泰勒展开得到：

f (r) + (Δ r) f^{'} (r) + ε g (r) + ε (Δ r) g^{'} (r) + O ({(Δ r)}^{2}) = 0

忽略高阶项，得到

(Δ r) (f^{'} (r) + ε g^{'} (r)) \approx - f (r) - ε g (r) = ε g (r)

或者，由于ε很小，所以εg'(r)很小。

Δ r \approx \frac{ε g (r)}{f^{'} (r) + ε g^{'} (r)} \approx \frac{ε g (r)}{f^{'} (r)}

根的敏感公式

假设r是函数f(x)的根，并且r+Δr是f(x)+εg(r)的根，则当ε<<f'(r)时，
$Δ r \approx \frac{ε g (r)}{f^{'} (r)}$

举例 $P(x) = (x - 1)(x - 2)(x - 3)(x - 4)(x - 5)(x - 6) - {10^{ - 6}}{x^7}$ 的最大根。

$f(x) = (x - 1)(x - 2)(x - 3)(x - 4)(x - 5)(x - 6)$ $\varepsilon = - {10^{ - 6}}$ $g(x)=x^7$ 。如果没有εg(x)，则有最大根r=6，但当加上这一项后，根变成了

r + Δ r = 6 - \frac{ε 6^{7}}{5!} = 6 - 2332.8 ε = 6.0023328

误差放大因子 = 相对前向误差 / 相对后向误差

条件数也是误差放大度量的一种方式，数值分析是对算法的研究，算法把定义问题的数据作为输入，对应的结果作为输出。条件数指的是理论问题本身所带来的误差放大部分，和用于求解问题的特定算法无关。注意到条件数仅仅度量由于问题本身带来的误差放大，这点很重要。和条件一起，还有一个平行概念，即稳定。稳定指的是由于算法小的输入误差造成的放大，而不是问题本身。如果一个算法在小的后向误差存在的时候，总能给出一个近似解，则称该算法是稳定的。如果问题的条件好，算法稳定，我们可以期望同时具有小的后向误差和前向误差。

前面的误差放大例子表明根求解过程对于特定的输入变化的敏感性。问题可能或多或少地敏感，依赖于如何设计输入的变化。问题的条件数定义为所有输入变化，或者至少规定类型的变化所造成的最大误差放大。条件数高的问题称为病态问题，条件数在1附近的问题称为良态问题。

牛顿方法

牛顿方法

$x_0$ =初始估计
$x_{i+1}=x_i-{f(x_i)\over f'(x_i)},\quad i=0,1,2,...$

定义 $e_i$ 表示一个迭代方法第i步后得到的误差，如果满足下式，则该迭代为二次迭代

M = lim_{i \to \infty} \frac{e_{i + 1}}{e_{i}^{2}} < \infty

定理 $e_i$ 满足

lim_{i \to \infty} \frac{e_{i + 1}}{e_{i}^{2}} = M

证明参见原书，注意此处

M = | \frac{f^{″} (r)}{2 f^{'} (r)} |

$x^2$ 的实根，此时牛顿迭代公式为

x_{i + 1} = \frac{x_{i}}{2}

$e_{i+1}=e_i/2$ ，是线性收敛。

定理 $e_i$ 满足

lim_{i \to \infty} \frac{e_{i + 1}}{e_{i}} = S = \frac{m - 1}{m}

定理：如果在[a,b]区间上f是(m+1)阶连续函数，包含m>1的多重根，则改进的牛顿方法为：

x_{i + 1} = x_{i} - \frac{m f (x_{i})}{f^{'} (x_{i})}

收敛到r，并具有二次收敛速度。

牛顿方法如同FPI（不动点迭代），可能不会收敛到根．下面的例子是一种可能不收敛的情况。

举例 $4x^4-6x^2-11/4$ $x_0=1/2$

不需要导数的根求解

除了重根，牛顿方法比二分法和FPI方法的收敛速度更快。它达到了这种更快的速度是因为使用了更多的信息，尤其是通过函数导数得到的函数切线方向的信息。在某些情况下，可能难以计算导数。在这种情况下，割线方法是牛顿方法的一个非常好的替代。它使用近似值割线替代了切线，并且收敛速度差不多快。

割线方法及其变体

$x_i$ 处的导数的近似可写为差商：

\frac{f (x_{i}) - f (x_{i - 1})}{x_{i} - x_{i - 1}}

割线方法

$x_0$ $x_1$ =初始估计
${x_{i + 1}} = {x_i} - {{f({x_i})\left( {{x_i} - {x_{i - 1}}} \right)} \over {f({x_i}) - f({x_{i - 1}})}},\quad i = 1,2,3, \cdots$

割线方法以超线性的速度收敛到一个单根，意味着它在线性和二次收敛方法之间。

试位方法和二分法相似，将中点从a和b的平均值改为割线方法定义的点，注意f(a)f(b)<0

c = a - \frac{f (a) (a - b)}{f (a) - f (b)} = \frac{b f (a) - a f (b)}{f (a) - f (b)}

c对应的便是a和b所定义的直线与x轴的交点。

试位方法开始表现得比二分法和割线方法都要好，具有二者最好的性质。但是，二分法在每一步中可以确保消除1/2的不确定性，试位方法却没有能力做出这样的保证，而且在一些例子中可能收敛很慢。

Muller方法 $x_0$ $x_1$ $x_2$ $x_2$ $x_3$ ，通过简单的二次公式计算，就可以确定两种可能。如果抛物线和x轴不相交，就会出现复数解，能够处理复数代数的软件就可以计算对应的解。

逆二次插值(IQI)是割线方法到抛物线的一种相近的泛化方法，但是使用形如x=p(y)的抛物线，而不是 Muller方法所使用的y=p(x)，这个抛物线和x轴只有一个交点。经过三点(a,A)，(b,B)，(c,C)的二阶多项式x=P(y)为：

P (y) = a \frac{(y - B) (y - C)}{(A - B) (A - C)} + b \frac{(y - A) (y - C)}{(B - A) (B - C)} + c \frac{(y - A) (y - B)}{(C - A) (C - B)}

其中P(A)=a，P(B)=b，P(C)=c，用y=0带入可得得到抛物线和x轴的交点，可以得到

P (0) = c - \frac{r (r - q) (c - b) + (1 - r) s (c - a)}{(q - 1) (r - 1) (s - 1)}

其中q=f(a)/f(b)，r=f(c)/f(b)，s=f(c)/f(a)。

逆二次插值

$a=x_i$ $b=x_{i+1}$ $c=x_{i+2}$ $A=f(x_i)$ $B=f(x_{i+1})$ $C=f(x_{i+2})$
$x_{i+3}=P(0)$ 为
${x_{i + 3}} = {x_{i + 2}} - {{r(r - q)({x_{i + 2}} - {x_{i + 1}}) + (1 - r)s({x_{i + 2}} - {x_i})} \over {(q - 1)(r - 1)(s - 1)}}$
$x_{i+3}$ $x_i$ ，另一种方式使用最新估计替换最近的三个估计中的后向误差最大的那个。

Brent方法

$x_i$ $[a_i,b_i]$ $x_i$ $a_i$ $b_i$ 中的一个：

（1）后向误差得到改进；

（2）包含根的区间至少减小一半；

否则，尝试使用割线方法以实现相同的目的。如果割线方法也失败了，则使用二分法，保证至少减少一半的不确定性。

方程组

高斯消去法不过多介绍

LU分解

可以将高斯消去法的消去步骤表示为LU分解，其中L是单位矩阵进行高斯消去操作后的下三角矩阵，U是原矩阵经过高斯消去得到的矩阵。

一旦知道L和U，问题Ax=b可以写作LUx=b，定义新的“辅助”向量c=Ux，则回代是个两步的过程

（1）对于方程Lc=b，求解c

（2）对于方程Ux=c，求解x

但并不是所有的矩阵都能进行LU分解，如

[\begin{matrix} 0 & 1 \\ 1 & 1 \end{matrix}]

无法被LU分解。

误差来源

在高斯消去法中有两个潜在的误差来源，第一个是病态问题的概念和解对于输入数据的敏感性相关，第二个是淹没，在大部分问题中可以通过一个简单的修正，称为部分主元，进行避免。

误差放大和条件数

定义 $x=(x_1,\cdots,x_n)$ ${\left\| x \right\|_\infty } = \max \left| {{x_i}} \right|$ ，i=1,...,n，即x所有元素的最大值。

定义 $x_a$ $x_a$ 后向误差 ${\left\| {b - A{x_a}} \right\|_\infty }$ 前向误差 ${\left\| {x - {x_a}} \right\|_\infty }$ 。

举例：找出近似解[-1,3.0001]的后向误差和前向误差，方程组如下：

\begin{aligned} x_{1} + x_{2} = 2 \\ 1.0001 x_{1} + x_{2} = 2.0001 \end{aligned}

精确解为[1,1]，后向误差计算：

b - A x_{a} = [\begin{matrix} 2 \\ 2.0001 \end{matrix}] - [\begin{matrix} 1 & 1 \\ 1.0001 & 1 \end{matrix}] [\begin{matrix} - 1 \\ 3.0001 \end{matrix}] = [\begin{matrix} - 0.0001 \\ 0.0001 \end{matrix}]

后向误差为0.0001，前向误差计算：

x - x_{a} = [\begin{matrix} 1 \\ 1 \end{matrix}] - [\begin{matrix} - 1 \\ 3.0001 \end{matrix}] = [\begin{matrix} 2 \\ - 2.0001 \end{matrix}]

前向误差为2.0001。可以看到前向误差和后向误差差别很大。下面这张图可以帮助理解这一情况

可以看到(1,1)是它们的交点，不过(-1,3.0001)也差一点在两条直线上（图中的差异是人为放大了1000倍，实际上离得非常近）。

$r=b-Ax_a$ ，系统Ax=b的相对后向误差定义为

\frac{{‖ r ‖}_{\infty}}{{‖ b ‖}_{\infty}}

相对前向误差定义为

\frac{{‖ x - x_{a} ‖}_{\infty}}{{‖ x ‖}_{\infty}}

误差放大因子则是二者的比值，即相对前向误差/相对后向误差。

定义：方阵A的条件数cond(A)为求解Ax=b时，对于所有右侧向量b，可能出现的最大误差放大因子。

方阵A的矩阵范数定义为每行元素绝对值求和，其中的最大值作为矩阵A的范数。

定理：n×n矩阵A的条件数是：

c o n d (A) = ‖ A ‖ \cdot ‖ A^{- 1} ‖

举例：计算矩阵A的条件数，其中

A = [\begin{matrix} 1 & 1 \\ 1.0001 & 1 \end{matrix}]

$\left\| A \right\|=2.0001$ ，A的逆为

A^{- 1} = [\begin{matrix} - 10000 & 10000 \\ 10001 & - 10000 \end{matrix}]

$\left\| A^{-1} \right\|=20001$ ，A的条件数是cond(A)=40004.0001，在这个系统中，对于任何其他的b，误差放大因子将小于或者等于40004.0001。

${\varepsilon _{mach}}$ ${\varepsilon _{mach}}$ ${\varepsilon _{mach}}$ $cond(A) \approx {10^k}$ ，我们在计算x时，将丢掉k位数字精度。

$cond(A) \approx 4 \times {10^4}$ ，对于双精度16位精确度而言，只能求解到大约12位的精确数字（实际上会更低，matlab只能求解到11位精确值）。

${H_{ij}} = 1/(i + j - 1)$ ），其对应的条件数非常大。

$\left\| {A + B} \right\| \le \left\| A \right\| + \left\| B \right\|$ 。

${\left\| x \right\|_1} = \left| {{x_1}} \right| + \cdots + \left| {{x_n}} \right|$ ${\left\| A \right\|_1}$ 的最大绝对列和，即列向量的1－范数的最大值。

淹没

举例：考虑方程组

\begin{aligned} 10^{- 20} x_{1} + x_{2} = 1 \\ x_{1} + 2 x_{2} = 4 \end{aligned}

（1）精确解约为[2,1]

（2）IEEE双精度，

[\begin{matrix} 10^{- 20} & 1 & | & 1 \\ 1 & 2 & | & 4 \end{matrix}] \overset{r_{2} - 10^{20} r_{1}}{⟶} [\begin{matrix} 10^{- 20} & 1 & | & 1 \\ 0 & 2 - 10^{20} & | & 4 - 10^{20} \end{matrix}]

${2 - {{10}^{20}}}$ ${4 - {{10}^{20}}}$ ${ - {{10}^{20}}}$ ，因此解为[0,1]

（3）IEEE双精度，使用行交换

[\begin{matrix} 1 & 2 & | & 4 \\ 10^{- 20} & 1 & | & 1 \end{matrix}] \overset{r_{2} - 10^{- 20} r_{1}}{⟶} [\begin{matrix} 1 & 2 & | & 4 \\ 0 & 1 - 2 \times 10^{- 20} & | & 1 - 4 \times 10^{- 20} \end{matrix}]

${1 - 2\times{{10}^{-20}}}$ ${1 - 4\times{{10}^{-20}}}$ 存为1，解得[2,1]，注意这里的结果是确切的结果，不像精确解存在微小的差距。

$10^{20}$ $10^{20}$ 倍，底部方程会被抑制，或者称为"swamp,"，尽管初始的时候有两个独立的方程或者源信息，但是经过第2种计算的消去，这里仅仅留下了顶部方程的两个副本。由于底部方程消失，出于所有实际的目的，我们不能指望计算结果可以满足底部的方程，实际上得到的解也不能满足对应的方程。

$10^{-20}$ ，在消去之后，原始的两个方程仍然存在，并且变为三角形式，结果是更加精确的近似解。

高斯消去的过程中保证乘子尽可能小，同时避免淹没。幸运的是，通过对朴素的高斯消去法的简单修正，可以使得高斯消去法中的乘子的绝对值不大于1。这种新的原则，包括在表格中明智的行交换，称为部分主元方法，在下一节中将进行详细的描述。

PA=LU分解

部分主元

$|a_{pk}|$ ，必要时交换第p行和第k行，然后继续进行消去。注意到使用部分主元方法保证所有乘子，或者L的元素的绝对值不大于1。通过这种对于高斯消去方法所进行的小的改变，淹没问题可以完全避免。

举例：使用高斯消去的部分主元方法求解方程组

\begin{aligned} x_{1} - x_{2} + 3 x_{3} = - 3 \\ - x_{1} - 2 x_{3} = 1 \\ 2 x_{1} + 2 x_{2} + 4 x_{3} = 0 \end{aligned}

$a_{31}$ 最大，所以交换第一行和第三行

\begin{aligned} [\begin{matrix} 1 & - 1 & 3 & | & - 3 \\ - 1 & 0 & - 2 & | & 1 \\ 2 & 2 & 4 & | & 0 \end{matrix}] \overset{r_{1} \leftrightarrow r_{3}}{⟶} [\begin{matrix} 2 & 2 & 4 & | & 0 \\ - 1 & 0 & - 2 & | & 1 \\ 1 & - 1 & 3 & | & - 3 \end{matrix}] \\ \overset{r_{2} + \frac{1}{2} r_{1}}{⟶} [\begin{matrix} 2 & 2 & 4 & | & 0 \\ 0 & 1 & 0 & | & 1 \\ 1 & - 1 & 3 & | & - 3 \end{matrix}] \overset{r_{3} - \frac{1}{2} r_{1}}{⟶} [\begin{matrix} 2 & 2 & 4 & | & 0 \\ 0 & 1 & 0 & | & 1 \\ 0 & - 2 & 1 & | & - 3 \end{matrix}] \end{aligned}

$a_{32}$ $a_{22}$ 大，所以交换第二行和第三行

\begin{aligned} [\begin{matrix} 2 & 2 & 4 & | & 0 \\ 0 & 1 & 0 & | & 1 \\ 0 & - 2 & 1 & | & - 3 \end{matrix}] \overset{r_{2} \leftrightarrow r_{3}}{⟶} [\begin{matrix} 2 & 2 & 4 & | & 0 \\ 0 & - 2 & 1 & | & - 3 \\ 0 & 1 & 0 & | & 1 \end{matrix}] \\ \overset{r_{3} + \frac{1}{2} r_{2}}{⟶} [\begin{matrix} 2 & 2 & 4 & | & 0 \\ 0 & - 2 & 1 & | & - 3 \\ 0 & 0 & 0.5 & | & - 0.5 \end{matrix}] \end{aligned}

解得x=[1,1,-1]

置换矩阵

定义：置换矩阵是一个n×n的矩阵，其在每一行、每一列仅有一个1，其他全部为0。

定理：令P是通过对单位矩阵实施一组特定的行交换后得到的一个n×n的置换矩阵。则对于任意的n×n矩阵 A，PA对应于对矩阵A实施同样的行交换得到的结果。

PA=LU分解

我们把关于高斯消去所知道的一切组合起来，进行PA=LU分解。这是部分主元的高斯消去的矩阵形式，PA=LU分解是求解线性方程组的主要方法。正如名字中所暗示的，PA=LU是对于矩阵A包含行交换的LU分解。在部分主元中，开始的时候我们并不知道需要进行置换的列，所以我们必须谨慎地把行置换的信息加入分解中，我们需要记录高斯消去法中所有的主元。

使用PA=LU求解方程组Ax=b时

\begin{aligned} P A x = P b \\ L U x = P b \\ L c = P b \Rightarrow c \\ U x = c \Rightarrow x \end{aligned}

其中P即为单位阵进行部分主元法交换行后产生的置换矩阵。

迭代方法

雅可比方法

雅可比(Jacobi)方法是方程组系统中的一种形式的不动点迭代。在FPI中，第一步是重写方程，进而求解未知量。雅可比方法以如下标准方式进行该重写步骤：求解第i个方程得到第i个未知变量，然后使用不动点迭代，从初始估计开始，进行迭代。

举例：使用雅可比方法求解下面的线性方程组

\begin{aligned} 3 x_{1} - x_{2} + 2 x_{3} = - 5 \\ - x_{1} + 4 x_{2} - 2 x_{3} = 4 \\ 2 x_{1} + x_{2} + 4 x_{3} = 0 \end{aligned}

求解未知变量得到

\begin{aligned} x_{1} = \frac{x_{2} - 2 x_{3} - 5}{3} \\ x_{2} = \frac{x_{1} + 2 x_{3} + 4}{4} \\ x_{3} = \frac{- 2 x_{1} - x_{2}}{4} \end{aligned}

迭代过程为：

\begin{aligned} [\begin{matrix} x_{1} \\ x_{2} \\ x_{3} \end{matrix}] = [\begin{matrix} 0 \\ 0 \\ 0 \end{matrix}] \\ [\begin{matrix} x_{1} \\ x_{2} \\ x_{3} \end{matrix}] = [\begin{matrix} \frac{x_{2} - 2 x_{3} - 5}{3} \\ \frac{x_{1} + 2 x_{3} + 4}{4} \\ \frac{- 2 x_{1} - x_{2}}{4} \end{matrix}] = [\begin{matrix} \frac{- 5}{3} \\ \frac{4}{4} \\ \frac{0}{4} \end{matrix}] = [\begin{matrix} \frac{- 5}{3} \\ 1 \\ 0 \end{matrix}] \\ [\begin{matrix} x_{1} \\ x_{2} \\ x_{3} \end{matrix}] = [\begin{matrix} \frac{x_{2} - 2 x_{3} - 5}{3} \\ \frac{x_{1} + 2 x_{3} + 4}{4} \\ \frac{- 2 x_{1} - x_{2}}{4} \end{matrix}] = [\begin{matrix} \frac{1 - 2 \times 0 - 5}{3} \\ \frac{\frac{- 5}{3} + 2 \times 0 + 4}{4} \\ \frac{- 2 \times \frac{- 5}{3} - 1}{4} \end{matrix}] = [\begin{matrix} \frac{- 4}{3} \\ \frac{7}{12} \\ \frac{7}{12} \end{matrix}] \\ . . . \\ [\begin{matrix} x_{1} \\ x_{2} \\ x_{3} \end{matrix}] = [\begin{matrix} - 1.833 \\ 0.889 \\ 0.694 \end{matrix}] \end{aligned}

定义 $A=(a_{ij})$ 是严格对角占优矩阵。

定理：如果n×n矩阵A是严格对角占优矩阵，则（1） A是非奇异矩阵，（2）对于所有向量b和初始估计，对Ax=b应用雅可比方法都会收敛到（唯一）解。

但是这个只是充分条件，不代表不满足严格对角占优矩阵，就不会收敛到唯一解，上例中就不是严格对角占优矩阵，可以看到第一行的对角线元素为3，只是等于第一行其它元素的绝对值之和。

雅可比方法是不动点迭代的一种形式，令D表示A的主对角线矩阵，L表示矩阵A的下三角矩阵（主对角线以下的元素），U表示上三角矩阵（主对角线以上的元素），则A=L+D+U，求解的方程变为Lx+Dx+Ux=b。注意，这里对于L和U的使用和LU分解中不同，当前L和U的主对角线元素都是零。

Ax=b

(D+L+U)x=b

Dx=b-(L+U)x

$D^{-1}$ (b-(L+U)x)

雅可比方法

$x_0$ =初始向量
$x_{k+1}$ $D^{-1}$ $x_k$ ), k = 0,1,2,...

高斯-塞德尔方法和SOR

高斯－塞德尔方法和雅可比方法的唯一差异是在前者中，最近更新的未知变量的值在每一步中都使用，即使更新发生在当前步骤。在使用相同步骤的情况下，该近似更加精确。如果收敛，高斯-塞德尔方法常常比雅可比方法收敛更快。

高斯-塞德尔方法

$x_0$ =初始向量
$x_{k+1}$ $D^{-1}$ $x_k$ $x_{k+1}$ ) k = 0,1,2,...

举例：对如下系统应用高斯-塞德尔方法

\begin{aligned} 3 u - v + 2 w = - 5 \\ - u + 4 v - 2 w = 4 \\ 2 u + v + 4 w = 0 \end{aligned}

迭代公式为

[\begin{matrix} u_{k + 1} \\ v_{k + 1} \\ w_{k + 1} \end{matrix}] = [\begin{matrix} \frac{v_{k} - 2 w_{k} - 5}{3} \\ \frac{u_{k + 1} + 2 w_{k} + 4}{4} \\ \frac{- 2 u_{k + 1} - v_{k + 1}}{4} \end{matrix}]

连续过松弛 $x_{k+1}$ $x_k$ 的平均。数字w被称为松弛参数，而当w>1时被称为过松弛。

使用SOR后的迭代公式

[\begin{matrix} u_{k + 1} \\ v_{k + 1} \\ w_{k + 1} \end{matrix}] = [\begin{matrix} (1 - w) u_{k} + w \frac{v_{k} - 2 w_{k} - 5}{3} \\ (1 - w) v_{k} + w \frac{u_{k + 1} + 2 w_{k} + 4}{4} \\ (1 - w) w_{k} + w \frac{- 2 u_{k + 1} - v_{k + 1}}{4} \end{matrix}]

可以将系统看成不动点迭代，如下：

\begin{aligned} A x = b \\ (w L + w D + w U) x = w b \\ (w L + D) x = w b - w U x + (1 - w) D x \\ x = {(w L + D)}^{- 1} [(1 - w) D x - w U x] + w {(w L + D)}^{- 1} b \end{aligned}

连续过松弛（SOR）

$x_0$ =初始向量
${x_{k + 1}} = {\left( {wL + D} \right)^{ - 1}}\left[ {\left( {1 - w} \right)D{x_k} - wU{x_k}} \right] + w{\left( {wL + D} \right)^{ - 1}}b$

收敛性的证明略

稀疏矩阵计算

$n^3$ $n^2$ 次的乘法（对于每个矩阵元素做一次乘法操作）以及大约相同数量的加法。

如果已知解的较好的近似，在这种特定的情况下支持使用迭代技术。例如，假设知道Ax＝b的解，随后A与b同时或者单个仅仅发生小的变化。我们可以想象一个动态系统，当A和b改变时，对A和b进行待续的度量，因而需要待续更新精确的解向量x。如果把前面问题的解作为新的相似问题的初始估计，使用雅可比或者高斯－塞德尔可以得到更快的收敛。这种技术通常称为修饰，这是由于该问题从一个近似解开始，该近似解可能对应前面一个相关问题的解，然后仅仅修饰近似解使其更加精确。

$n^2$ 个矩阵元素，只有O(n)个非零元素。一个完全的矩阵恰恰相反，其中没有几个元素可以假设为0。对稀疏矩阵使用高斯消去经常会导致填充，这使得其中系数矩阵由于必要的行交换从稀疏变为不稀疏。由于这个原因，高斯消去以及PA=LU实现的效率对于稀疏矩阵变得可疑，在这种情况下迭代方法是一种合理的选择。


x
% 程序 2.1 稀疏矩阵生成
% 输入：n = 系统大小
% 输出：稀疏矩阵 a, r.h.s. b
function [a,b] = sparsesetup (n)
e = ones(n,l); n2=n/2;
a= spdiags( [-e 3*e -e], -1 : 1,n,n);   % a的元素
c= spdiags( e/2, 0,n,n) ; c= fliplr(c); a= a+c;
a(n2+1,n2) = -1; a (n2,n2+1) = -1;      % 设置2个元素
b=zeros(n,1);                           % r.h.s. b的元素
b(1)=2.5; b(n)=2.5; b(2:n-1)=1.5; b(n2:n2+1)=1;

% 程序 2.2 雅可比方法
% 输入：完全矩阵或者稀疏矩阵 a, r.h.s. b,
% 雅可比迭代次数k
% 输出：解x
function x = jacobi(a, b, k)
n=length(b);            % 确定n
d=diag(a);              % 提取a的对角线元素
r=a-diag(d );           % r为余项
x=zeros(n,1);           % 初始化向量x
for j =1:k              % 雅可比迭代的循环
    x = (b-r*x) ./ d;
end

用来对称正定矩阵的方法

对称正定矩阵

定义 $A^T=A$ 对称矩阵 $x^TAx>0$ ，则矩阵A是正定矩阵。

性质：如果n×n矩阵A是对称矩阵，则A是正定矩阵，当且仅当所有特征值是正数。

性质 $X^TAX$ 是m×m对称正定矩阵。

定义：方阵A的主子矩阵是一个方的子矩阵，其对角线元素是矩阵A的对角线元素。

性质：任何对称正定矩阵的主子矩阵也是对称正定矩阵。

Cholesky分解

考虑一个对称正定矩阵

[\begin{matrix} a & b \\ b & c \end{matrix}]

$R^TR$ 的形式

[\begin{matrix} a & b \\ b & c \end{matrix}] = [\begin{matrix} \sqrt{a} & 0 \\ u & v \end{matrix}] [\begin{matrix} \sqrt{a} & u \\ 0 & v \end{matrix}] = [\begin{matrix} a & u \sqrt{a} \\ u \sqrt{a} & u^{2} + v^{2} \end{matrix}]

$u = b/\sqrt a$ $v^2=c-u^2$ ，Cholesky分解如下：

A = [\begin{matrix} a & b \\ b & c \end{matrix}] = [\begin{matrix} \sqrt{a} & 0 \\ \frac{b}{\sqrt{a}} & \sqrt{c - b^{2} / a} \end{matrix}] [\begin{matrix} \sqrt{a} & \frac{b}{\sqrt{a}} \\ 0 & \sqrt{c - b^{2} / a} \end{matrix}] = R^{T} R

定理 $A=R^TR$ 。

Cholesky分解

for k = 1,2,...,n
$A_{kk}$ <0，stop，end
$R_{kk}=\sqrt{A_{kk}}$
$u^T={1\over R_{kk}}A_{k,k+1:n}$
$R_{k,k+1:n}$ $u^T$
$A_{k+1:n,k+1:n}=A_{k+1:n,k+1:n}-uu^T$
end

举例：计算如下矩阵的Cholesky分解

[\begin{matrix} 4 & - 2 & 2 \\ - 2 & 2 & - 4 \\ 2 & - 4 & 11 \end{matrix}]

$R_{11}=\sqrt{A_{11}}=2$ $u^T=A_{1,2:3}/R_{11}=[-2,2]/2=[-1,1]$

$R_{1,2:3}=[-1,1]$ $A_{2:3,2:3}$ 进行运算：

[\begin{matrix} 2 & - 4 \\ - 4 & 11 \end{matrix}] - [\begin{matrix} 1 & - 1 \\ - 1 & 1 \end{matrix}] = [\begin{matrix} 1 & - 3 \\ - 3 & 10 \end{matrix}]

$R_{22}=\sqrt{A_{22}}=1$ $u^T=A_{23}=-3$ $R_{23}=u^T=-3$ $R_{33}=\sqrt{1}$ ，所以R等于

R = [\begin{matrix} 2 & - 1 & 1 \\ 0 & 1 & - 3 \\ 0 & 0 & 1 \end{matrix}]

共轭梯度方法

$(v,w)=v^Tw$ 对称并对于输入v和w线性，欧几里得内积也是正定的，当v≠0时，(v,v)>0。

定义：令A是对称正定的n×n矩阵，对于两个n维向量v和w，定义A内积

(v, w)_{A} = v^{T} A w

$(v,w)_A=0$ 时，向量v和w为A共轭。

共轭梯度方法

$x_0$ =初始估计
$d_0=r_0=b-Ax_0$
for k = 0,1,2,...,n-1
$r_k$ =0，stop，end
${\alpha _k} = {{{r_k}^T{r_k}} \over {d_k^TA{d_k}}}$
$x_{k+1}=x_k+\alpha_kd_k$
$r_{k+1}=r_k-\alpha_kAd_k$
$\beta_k={r_{k+1}^Tr_{k+1}\over r_k^Tr_k}$
$d_{k+1}=r_{k+1}+\beta_kd_k$
end

$x_k$ $r_k$ $d_k$ $x_k$ $x_{k+1}$ $d_k$ ，并使之两两共扼。共轭的概念推广了正交的概念，并据此在算法的名字中也包含“共轭”。

$\alpha_k$ $\beta_k$ $\alpha_k$ $r_{k+1}$ $d_k$ 正交：

\begin{aligned} d_{k}^{T} r_{k + 1} = d_{k}^{T} r_{k} - α_{k} d_{k}^{T} A d_{k} = 0 \\ α_{k} = \frac{{r_{k}}^{T} r_{k}}{d_{k}^{T} A d_{k}} \end{aligned}

$\beta_k$ $d_{k-1}$ $r_k$ 正交，有

\begin{aligned} d_{k} - r_{k} = β_{k - 1} d_{k - 1} \\ r_{k}^{T} d_{k} - r_{k}^{T} r_{k} = 0 \end{aligned}

$r_k^T{d_k} = r_k^T{r_k}$ $d_k$ 两两A共轭，

\begin{aligned} d_{k + 1} = r_{k + 1} + β_{k} d_{k} \\ d_{k}^{T} A d_{k + 1} = d_{k}^{T} A r_{k + 1} + β_{k} d_{k}^{T} A d_{k} \\ 0 = d_{k}^{T} A r_{k + 1} + β_{k} d_{k}^{T} A d_{k} \\ β_{k} = - \frac{d_{k}^{T} A r_{k + 1}}{d_{k}^{T} A d_{k}} \end{aligned}

$\beta_k$ 还可以继续化简成伪代码中的形式，具体的证明参考原书。

预条件

通过使用预条件技术，可以使得诸如共扼梯度方法的迭代方法收敛速度加快。迭代方法的收敛率通常直接或者间接依赖于系数矩阵A的条件数。预条件方法的思想是降低问题中的条件数。

n×n的线性方程组Ax=b的预条件形式是

M^{- 1} A x = M^{- 1} b

其中M是可逆的n×n矩阵，称为预条件子。我们所要做的就是在方程两侧左乘上该矩阵，预条件矩阵试图对矩阵A逆转，从而可以有效地降低问题的条件数。概念上来讲，它试图同时做两件事：矩阵M应该（1）和矩阵A足够接近，（2）容易求逆。这两个目的常常彼此对立。

一种特别简单的方式是雅可比预条件子M=D，其中D是A的对角线矩阵。D的逆矩阵是D的元素的倒数．例如在一个严格对角占优矩阵中，雅可比预条件子和A相似，同时非常容易求逆。注意到，对称正定矩阵的每个对角线元素都严格为正，因而计算倒数不是问题。

预条件共轭梯度法

$x_0$ =初始估计
$r_0=b-A_0$
$d_0=z_0=M^{-1}r_0$
for k = 0,1,2,...,n-1
$r_k$ =0，stop，end
$a_k=r_k^Tz_k/d^T_kAd_k$
$x_{k+1}=x_k+\alpha_kd_k$
$r_{k+1}=r_k-\alpha_kAd_k$
$z_{k+1}=M^{-1}r_{k+1}$
$\beta_k=r_{k+1}^Tz_{k+1}/r_k^Tz_k$
$d_{k+1}=z_{k+1}+\beta_kd_k$
end

在对称连续过松弛(SSOR)中，预条件子定义如下：

M = (D + w L) D^{- 1} (D + w U)

其中A=L+D+U被分为下三角部分、对角线以及上三角部分。w是一个在0和2之间的常数。在特例中w=1，这被称为高斯-塞德尔预条件子。

$M=(I+wLD^{-1})(D+wU)$ $z=M^{-1}v$ 可以通过两次回代求解。

\begin{matrix} (I + w L D^{- 1}) c = v \\ (D + w U) z = c \end{matrix}

非线性方程组

多元牛顿方法

多变量情况下函数导数f'对应的是雅可比矩阵，定义为

D F (x) = [\begin{matrix} \frac{\partial f_{1}}{\partial u} & \frac{\partial f_{1}}{\partial v} & \frac{\partial f_{1}}{\partial w} \\ \frac{\partial f_{2}}{\partial u} & \frac{\partial f_{2}}{\partial v} & \frac{\partial f_{2}}{\partial w} \\ \frac{\partial f_{3}}{\partial u} & \frac{\partial f_{3}}{\partial v} & \frac{\partial f_{3}}{\partial w} \end{matrix}]

多变量牛顿方法

$x_0$ =初始向量
$x_{k+1}=x_k-(DF(x_k))^{-1}F(x_k)$ ，k=0,1,2,...

$x_{k+1}=x_k-s$ $DF(x_k)s=F(x_k)$ 的解，所以迭代步骤为

{\begin{matrix} D F (x_{k}) s = - F (x_{k}) \\ x_{k + 1} = x_{k} + s \end{matrix}

Broyden方法

如果可以计算雅可比矩阵，那么牛顿方法是一个好选择，那么最好的替代方法就是Broyden方法。

Broyden方法

$x_0$ =初始向量
$A_0$ =初始矩阵
for i = 0,1,2,...
$x_{i+1}=x_i-A_i^-F(x_i)$
${A_{i + 1}} = {A_i} + {{\left( {{\Delta _{i + 1}} - {A_i}{\delta _{i + 1}}} \right){\delta _{i + 1}}^T} \over {{\delta _{i + 1}}^T{\delta _{i + 1}}}}$
end

${\delta _{i + 1}} = {x_{i + 1}} - {x_i},\;{\Delta _{i + 1}} = F({x_{i + 1}}) - F({x_i})$

Broyden方法Ⅱ

$x_0$ =初始向量
$A_0$ =初始矩阵
for i = 0,1,2,...
$x_{i+1}=x_i-B_iF(x_i)$
${B_{i + 1}} = {B_i} + {{\left( {{\delta _{i + 1}} - {B_i}{\Delta _{i + 1}}} \right){\delta _{i + 1}}^T{B_i}} \over {{\delta _{i + 1}}^T{B_i}{\delta _{i + 1}}}}$
end

${\delta _{i + 1}} = {x_{i + 1}} - {x_i},\;{\Delta _{i + 1}} = F({x_{i + 1}}) - F({x_i})$

$x_0$ $B_0$ $B_0=I$ $B_i$ $A_i$ 估计雅可比。两个版本的Broyden方法都超线性收敛到单根，比牛顿方法的二次收敛要慢一些。

插值

数据和插值函数

定义 $P(x_i)=y_i$ $(x_1,y_1)$ $(x_n,y_n)$ 。

拉格朗日插值

$x_1$ $y_1$ $x_2$ $y_2$ $x_3$ $y_3$ )，则多项式

P_{2} (x) = y_{1} \frac{(x - x_{2}) (x - x_{3})}{(x_{1} - x_{2}) (x_{1} - x_{3})} + y_{2} \frac{(x - x_{1}) (x - x_{3})}{(x_{2} - x_{1}) (x_{2} - x_{3})} + y_{3} \frac{(x - x_{1}) (x - x_{2})}{(x_{3} - x_{1}) (x_{3} - x_{2})}

$(x_1,y_1)$ $(x_n,y_n)$ ，对于1~n之间的每个k，定义n-1次多项式

L_{k} (x) = \frac{(x - x_{1}) \dots (x - x_{k - 1}) (x - x_{k + 1}) \dots (x - x_{n})}{(x_{k} - x_{1}) \dots (x_{k} - x_{k - 1}) (x_{k} - x_{k + 1}) \dots (x_{k} - x_{n})}

${L_k}(x)$ ${L_k}(x_k)=1$ ${L_k}(x_j)=0\;j\neq k$ 。定义了n-1次多项式

P_{n - 1} (x) = y_{1} L_{1} (x) + \dots y_{n} L_{n} (x)

定理（多项式插值的主定理） $(x_1,y_1)$ $(x_n,y_n)$ $x_i$ $P(x_i)=y_i$ ，i=1,...,n。

牛顿差商

定义 $x_1$ $x_n$ $x^{n-1}$ $x_1$ $x_1$ $x_n$ $x_n$ ))。

牛顿差商公式：

\begin{aligned} P (x) = f [x_{1}] + f [x_{1} x_{2}] (x - x_{1}) + f [x_{1} x_{2} x_{3}] (x - x_{1}) (x - x_{2}) \\ + f [x_{1} x_{2} x_{3} x_{4}] (x - x_{1}) (x - x_{2}) (x - x_{3}) + . . . \\ + f [x_{1} \dots x_{4}] (x - x_{1}) \dots (x - x_{n - 1}) \end{aligned}

牛顿差商

$x_1$ $x_n$ $y_1$ $y_n$ ]
for j = 1,...,n
$x_j$ $y_j$
end
for i = 2,...,n
for j=1,...,n+1-i
$x_j$ $x_{j+i-1}$ $x_{j+1}$ $x_{j+i-1}$ $x_{j}$ $x_{j+i-2}$ $x_{j+i-1}$ $x_j$ )
end
end

插值代码


xxxxxxxxxx
% 程序 3.1 牛顿差商插值方法
% 计算插值多项式的系数
% 输入：x 和y是包含 n 个数据点的 x 和 y坐 标的向量
% 输出：嵌套形式的插值多项式系数 c
% 使用 nes 七 ． m 计算插值 多项式
function c=newtdd(x, y,n)
for j =1:n
    v(j,1) =y (j) ;     % 填入牛顿三角形的y列
end
for i =2:n  % 对于第l.列
    for j =1 : n+1- i % 从顶端到底端填入列元素
        v (j, i) = (v (j + 1, i -1) - v(j, i -1)) / (x (j +i - 1) -x (j));
    end
end
for i = 1 : n
    c(i) =v(1, i) ; % 从三角形顶端读
end  % 输出系数

插值误差

不过多赘述。

切比雪夫插值

切比雪夫插值是一种特定最优的点间距选取方式。

切比雪夫理论

切比雪夫插值的动机是在插值区间上，提高对如下插值误差的最大值的控制

\frac{(x - x_{1}) (x - x_{2}) \dots (x - x_{n})}{n!} f^{(n)} (c)

从现在开始让我们把区间固定在[-1,1]。多项式插值误差的分子

(x - x_{1}) (x - x_{2}) \dots (x - x_{n})

$x_1$ $x_n$ 使得多项式插值误差的分子的最大值足够小？这被称为插值误差的最小最大问题。

定理 $-1\le x_1,\cdots, x_n\le 1$ ，使得

max_{- 1 \leq x \leq 1} | (x - x_{1}) \dots (x - x_{n}) |

尽可能小，则

x_{i} = \cos \frac{(2 i - 1) π}{2 n}, i = 1, \dots n

$1/2^{n-1}$ ，实际上，通过

(x - x_{1}) \dots (x - x_{n}) = \frac{1}{2^{n - 1}} T_{n} (x)

$T_n(x)$ 表示n阶切比雪夫多项式。

切比雪夫多项式

$T_n(x)=\cos (n\arccos x)$ ，对于每个n都是关于x的多项式。令y=arccosx，因而cosy=x。

\begin{aligned} T_{n + 1} (x) = \cos (n + 1) y = \cos (n y + y) = \cos n y \cos y - \sin n y \sin y \\ T_{n - 1} (x) = \cos (n - 1) y = \cos (n y - y) = \cos n y \cos y + \sin n y \sin y \\ T_{n + 1} (x) + T_{n - 1} (x) = 2 \cos n y \cos y = 2 x T_{n} (x) \\ T_{n + 1} (x) = 2 x T_{n} (x) - T_{n - 1} (x) \end{aligned}

事实 $\deg ({T_n}) = n$ $2^{n-1}$ 。

事实 $T_n(1)=1$ $T_n(-1)=(-1)^n$ 。

事实 $T_n(x)$ $T_n(x)=\cos y$ 的形式得到。

事实 $T_n(x)$ 的所有过零点都在-1~1之间，实际上，过零点是0=cos(n arccosx)的解，由于cosy=0，当且仅当y=奇数∙（π／2）。

\begin{aligned} n \arccos x = o d d \cdot π / 2 \\ x = \cos \frac{o d d \cdot π}{2 n} \end{aligned}

事实 $T_n(x)$ 在-1和1之间一共往返变化n+1次，实际上，这发生在cos0，cosπ／n，...，cos(n-1)π/n，cosπ。

区间的变化

关于切比雪夫插值的讨论局限在区间[-1,1]，可以将方法推广到一般的区间[a,b]，移动基点使得它们在区间[a,b]上的相对位置和在区间[-1,1]上一致。这可以通过如下两步实现：

（1）使用因子(b-a)/2拉伸点（这是两个区间长度的比值），

（2）将点平移(b+a)／2，使得中心从0移动到区间[a,b]的中心。

换句话讲，从原始点

\cos \frac{o d d π}{2 n} \to \frac{b - a}{2} \cos \frac{o d d π}{2 n} + \frac{b - a}{2}

$x_1$ $x_n$ $x-x_i$ $2^{n-1}$ ${\left[ {\left( {b - a} \right)/2} \right]^n}/{2^{n - 1}}$ 。

切比雪夫插值节点

在区间[a,b]
$x_{i} = \frac{b + a}{2} + \frac{b - a}{2} \cos \frac{(2 i - 1) π}{2 n}$
i=1,...,n，不等式
$| (x - x_{1}) \dots (x - x_{n}) | \leq \frac{{(\frac{b - a}{2})}^{n}}{2^{n - 1}}$
在区间[a,b]上成立。

举例：在区间[0,π/2]找出4个切比雪夫基点进行插值，找出切比雪夫插值误差的上界，在区间中f(x)=sinx。

x_{i} = \frac{\frac{π}{2} + 0}{2} + \frac{\frac{π}{2} - 0}{2} \cos \frac{(2 i - 1) π}{8}, i = 1, 2, 3, 4

\begin{aligned} x_{1} = \frac{π}{4} + \frac{π}{4} \cos \frac{π}{8} x_{2} = \frac{π}{4} + \frac{π}{4} \cos \frac{3 π}{8} \\ x_{3} = \frac{π}{4} + \frac{π}{4} \cos \frac{5 π}{8} x_{4} = \frac{π}{4} + \frac{π}{4} \cos \frac{7 π}{8} \end{aligned}

插值误差的最坏情况是

\begin{aligned} | \sin x - P_{3} (x) | = \frac{| (x - x_{1}) (x - x_{2}) (x - x_{3}) (x - x_{4}) |}{4!} | f^{⁗} (c) | \\ \leq \frac{{(\frac{\frac{π}{2} - 0}{2})}^{4}}{4! 2^{3}} 1 \approx 0.00198 \end{aligned}

牛顿差商 $x_1$ $x_2$ ，...。

三次样条

样条是另一种数据插值的方式，在多项式插值中，多项式给出的单一公式满足所有数据点，而样条则使用多个公式，其中每个都是低阶多项式，来通过所有数据点。线性样条可以成功地对任意的n 个点集进行插值，但是线性样条函数缺乏平滑，三次样条则可以解决线性样条的这个缺点。三次样条在两个数据点之间使用3阶(cubic)多项式替换线性样条两点之间的线性函数。由于通过n个点的三次样条无穷多，所以需要添加额外条件。

样条的性质

$(x_1,y_1)$ $(x_n,y_n)$ $x_i$ 不同且升序。通过n个点的三次样条S(x)是一组三次多项式

\begin{aligned} S_{1} (x) = y_{1} + b_{1} (x - x_{1}) + c_{1} (x - x_{1})^{2} + d_{1} (x - x_{1})^{3} x \in [x_{1}, x_{2}] \\ S_{2} (x) = y_{2} + b_{2} (x - x_{2}) + c_{2} (x - x_{2})^{2} + d_{2} (x - x_{2})^{3} x \in [x_{2}, x_{3}] \\ ⋮ \\ S_{n - 1} (x) = y_{n - 1} + b_{n - 1} (x - x_{n - 1}) + c_{n - 1} (x - x_{n - 1})^{2} + d_{n - 1} (x - x_{n - 1})^{3} x \in [x_{n - 1}, x_{n}] \end{aligned}

并具有如下性质：

$S_i(x_i)=y_i, S_i(x_{i+1})=y_{i+1}$ ，i=1,2,...,n-1。

$S'_{i-1}(x_i)=S'_i(x_i)$ ，其中i=2,...,n-1。

$S''_{i-1}(x_i)=S''_i(x_i)$ ，其中i=2,...,n-1。

性质1保证样条S(x)插值数据点；性质2使得相邻的样条段在它们相遇的地方斜率相同；性质3则保证在两条样条段相邻的地方曲率相同，该曲率由二阶导数表示。

只有满足上面的这些条件才被称为三次样条

性质4a（自然样条） $S''_{1}(x_1)=S''_{n-1}(x_n)=0$

满足性质4a这两个附加条件的三次样条被称为自然三次样条，自然三次样条是唯一的。

$\delta_i=x_{i+1}-x_i$ $\Delta_i=y_{i+1}-y_i$

自然三次样条

$x=[x_1,\cdots,x_n]$ $x_1<\cdots<x_n$ $y=[y_1,\cdots,y_n]$
for i = 1,...,n-1
$a_i=y_i$
$\delta_i=x_{i+1}-x_i$
$\Delta_i=y_{i+1}-y_i$
end
$c_1,\cdots,c_n$
for i = 1,...,n-1
$d_i={c_{i+1}-c_i\over 3\delta_i}$
$b_i={\Delta_i\over \delta_i}-{\delta_i\over 3}(2c_{i}+c_{i+1})$
end

除了性质4a给出的自然三次样条，还有曲率调整三次样条、钳制三次样条、抛物线端点的三次样条、非纽结三次样条。

贝塞尔曲线

贝塞尔样条是一个允许用户控制节点处斜率的样条。作为额外自由控制的代价，在节点处的一阶和二阶导数的平滑性不再能保证，而这种平滑性是前面三次样条本身就具有的性质。贝塞尔样条适合不时出现角点（一阶导数不连续）和曲率突变（二阶导数不连续）的情况。

$(x_1,y_1)$ $(x_2,y_2)$ $(x_3,y_3)$ $(x_4,y_4)$ $(x_2-x_1,y_2-y_1)$ $(x_1,y_1)$ $(x_4-x_3,y_4-y_3)$ $(x_4,y_4)$ 点结束。

贝塞尔曲线

$(x_1,y_1)$ $(x_4,y_4)$
$(x_2,y_2)$ $(x_3,y_3)$
设
$\begin{aligned} b_{x} = 3 (x_{2} - x_{1}) \\ c_{x} = 3 (x_{3} - x_{2}) - b_{x} \\ d_{x} = x_{4} - x_{1} - b_{x} - c_{x} \\ b_{y} = 3 (y_{2} - y_{1}) \\ c_{y} = 3 (y_{3} - y_{2}) - b_{y} \\ d_{y} = y_{4} - y_{1} - b_{y} - c_{y} \end{aligned}$
定义在0≤t≤1的贝塞尔曲线如下：
$\begin{aligned} x (t) = x_{1} + b_{x} t + c_{x} t^{2} + d_{x} t^{3} \\ y (t) = y_{1} + b_{y} t + c_{y} t^{2} + d_{y} t^{3} \end{aligned}$

$x_1$ $y_1$ $x_2$ $y_2$ $x_3$ $y_3$ $x_4$ $y_4$ 来定义每一段的贝塞尔样条。

237 620 237 620 237 120 237 120; 237 120 237 35 226 24 143 19; 143 19 143 19 143 0 143 0; 143 0 143 0 435 0 435 0; 435 0 435 0 435 19 435 19; 435 19 353 23 339 36 339 109; 339 109 339 108 339 620 339 620; 339 620 339 620 393 620 393 620; 393 620 507 620 529 602 552 492; 552 492 552 492 576 492 576 492; 576 492 576 492 570 662 570 662; 570 662 570 662 6 662 6 662; 6 662 6 662 0 492 0 492; 0 492 0 492 24 492 24 492; 24 492 48 602 71 620 183 620; 183 620 183 620 237 620 237 620;

画出来的曲线为：

python代码： BezierCurve[贝塞尔曲线].py

最小二乘

最小二乘与法线方程

在方程组一章中，给出了当方程解存在时，如何找到Ax＝b的解。我们将学习当解不存在的时候该怎么办。当方程不一致时，有可能方程的个数超过未知变量的个数，答案是找到第二可能好的解：即最小二乘近似．

在插值一章中给出了如何找出多项式，并精确拟合数据点。但是如果有大量的数据点，或者采集的数据点具有一定误差，使用高阶多项式精确拟合一般不是个好方法。在这种情况下，使用简单模型近似数据是一种更合理的方式。

求解不一致的方程组以及近似拟合数据这两个问题一同驱动着最小二乘的发展．

不一致的方程组

如果一个方程组无解，它被称为不一致。方程无解意味着什么？可能系数不十分精确，在很多情况下，方程的个数比未知变量的个数要多，使得解不可能满足所有的方程。但是还可以通过找出与解最相似的向量x。

如果我们选择“相似度”意味着欧氏距离相近，有一个直接的方法找到最接近的x，这个特殊的x称为最小二乘解。

对于Ax=b这个方程而言，如果b不在Ax可以确定的平面（A中的每一列看成基向量，对应的x的一行（就是一个值）便可看成该方向的值）上，那么就是无解的，不过可以在这个平面上找到与b最接近的向量来近似。

最小二乘基于正交，从一点到一个平面的最短距离，由一个到平面的正交线段表示，法线方程可以确定该线段，这表示最小二乘的误差。

$\bar x$ $b-A\bar x$ 代表余弦误差，让误差与平面正交，即

\begin{aligned} A^{T} (b - A \bar{x}) = 0 \\ A^{T} A \bar{x} = A^{T} b \end{aligned}

此时这个方程便是可解的，不过解出的结果与实际结果可能存在一定差距，有三种方式可以表示余项的大小。

$||r||_2=\sqrt{r_1^2+\cdots+r_m^2}$

$SE=r_1^2+\cdots+r_m^2$

$RMSE=\sqrt{SE/m}=\sqrt{(r_1^2+\cdots+r_m^2)/m}$

三种表达紧密相关。

数据的拟合模型

最小二乘数据拟合

$(t_1,y_1)$ $(t_m,y_m)$
step1 选择模型，确定用于拟合数据的参数模型，例如y=a+bt
$c_1$ $c_2$ 为参数，这得到系统Ax=b，其中未知变量x表示未知参数。
${A^T}A\bar x = {A^T}b$ 的解。

最小二乘是数据压缩的经典例子，输入包含一组数据点，输出是一个尽可能好的数据拟合模型，其中具有较少的参数。通常，使用最小二乘的原因是使用合理的底层模型替换噪声数据，该模型通常用于信号预测以及分类。

最小二乘的条件

$A^TA$ 的条件数过大，难以在双精度算术中处理，即使原始的最小二乘问题条件还可以，但是这个法线方程是病态问题。

模型概述

周期数据

周期数据使用周期模型。例如外层大气温度遵循大时间尺度的循环，包含由地球自转和绕太阳公转控制的每天和每年的循环。

$y=c_1+c_2\cos 2\pi t+c_3 \sin 2\pi t + c_4 \cos 4\pi t$ 。

数据线性化

以描述人口变化的指数模型为例：

y = c_{1} e^{c_{2} t}

$c_2$ 在模型方程中不是线性的。一旦将数据点代入模型，求解的困难很明显：要求解的方程系数不是线性，不能表示为线性方程组Ax＝b。因而我们导出的法线方程不能使用。

不过有两种方法可以用来解决这个问题，第一种方法是直接求解非线性最小二乘方程（在后面介绍），第二种是通过取对数来使问题线性化。

\ln y = \ln (c_{1} e^{c_{2} t}) = \ln c_{1} + c_{2} t

$y=c_1e^{c_2t}$ ，也可以采用这种方法解决。

QR分解

格拉姆-施密特正交与最小二乘

格拉姆－施密特方法是对一组向量正交化。给定一组输入的m维向量，目的是找出正交坐标系统，获取由这些向量张成的空间。更精确地讲，给定n个线性无关的输入向量，该方法计算n个彼此垂直的单位向量（单位长度由2范数定义），张成和输入向量相同的子空间。

$r_{ij}=q_i^TA_j$ $q_i$ $r_{jj}=||y_j||_2$ $r_{ij}=q_i^TA_j$ ，则可以将格拉姆-施密特正交换一种方式表达

\begin{aligned} A_{1} = r_{11} q_{1} \\ A_{2} = r_{12} q_{1} + r_{22} q_{2} \\ \dots \\ A_{j} = r_{1 j} q_{1} + \dots + r_{j - 1, j} q_{j - 1} + r_{j j} q_{j} \end{aligned}

或者写成矩阵形式

(\begin{matrix} A_{1} & \dots & A_{n} \end{matrix}) = (\begin{matrix} q_{1} & \dots & q_{n} \end{matrix}) [\begin{matrix} r_{11} & r_{12} & \dots & r_{1 n} \\ r_{22} & \dots & r_{2 n} \\ ⋱ & ⋮ \\ r_{n n} \end{matrix}]

$A_j$ 消减QR分解 $A_j$ $r_{jj}$ 非0。

经典格拉姆-施密特正交

$A_j$ （j=1，...，n）为线性无关向量。
for j = 1，2，...，n
$y=A_j$
for i = 1，2，...，j-1
$r_{ij}=q_i^TA_j$
$y=y-r_{ij}q_i$
end
$r_{jj}=||y||_2$
$q_j=y/r_{jj}$
end

$R^m$ $A_j$ $R^m$ ，并实现格拉姆-施密特方法。如下所示

(\begin{matrix} A_{1} & \dots & A_{n} \end{matrix}) = (\begin{matrix} q_{1} & \dots & q_{m} \end{matrix}) [\begin{matrix} r_{11} & r_{12} & \dots & r_{1 n} \\ r_{22} & \dots & r_{2 n} \\ ⋱ & ⋮ \\ \begin{matrix} 0 \\ ⋮ \\ 0 \end{matrix} & \begin{matrix} \dots \\ \dots \end{matrix} & \begin{matrix} \dots \\ \dots \end{matrix} & \begin{matrix} r_{n n} \\ 0 \\ ⋮ \\ 0 \end{matrix} \end{matrix}]

注意A是m×n矩阵，Q是m×m方阵，上三角矩阵R是m×n矩阵，m>n，这称为A的完全QR分解，Q在数值分析中具有特殊的地位，并具有一个特殊的定义。

定义 $Q^T=Q^{-1}$ ，方阵Q正交。

引理：如果Q是m×m正交矩阵，x是m维向量，则

| | Q x | |_{2} = | | x | |_{2}

LU分解是对高斯消去中信息进行有效编码的方式。以相同方式，QR分解记录了矩阵正交化的信息，即构造一个正交集，张出由A的列向量构成的空间。使用正交矩阵计算更好的原因是

（1）根据定义它们很容易求逆，

（2）由引理知，它们不会放大误差。

QR分解可用于求解n个方程n个未知量的方程组，但是计算复杂度比LU分解大。

$||Ax-b||_2$ ，使用上面的引理得到：

\begin{aligned} {‖ Q x ‖}_{2} = {‖ x ‖}_{2} \\ {‖ Q R x - b ‖}_{2} = {‖ Q^{- 1} (Q R x - b) ‖}_{2} = {‖ R x - Q^{T} b ‖}_{2} \end{aligned}

欧几里得范数中的向量是

[\begin{matrix} e_{1} \\ ⋮ \\ e_{n} \\ e_{n + 1} \\ ⋮ \\ e_{m} \end{matrix}] = [\begin{matrix} r_{11} & r_{12} & \dots & r_{1 n} \\ r_{22} & \dots & r_{2 n} \\ ⋱ & ⋮ \\ \begin{matrix} 0 \\ ⋮ \\ 0 \end{matrix} & \begin{matrix} \dots \\ \dots \end{matrix} & \begin{matrix} \dots \\ \dots \end{matrix} & \begin{matrix} r_{n n} \\ 0 \\ ⋮ \\ 0 \end{matrix} \end{matrix}] [\begin{matrix} x_{1} \\ ⋮ \\ x_{n} \end{matrix}] - [\begin{matrix} d_{1} \\ ⋮ \\ d_{n} \\ d_{n + 1} \\ ⋮ \\ d_{m} \end{matrix}]

$d=Q^Tb$ $r_{ii}\neq0$ $e$ $(e_1,\cdots,e_n)$ $x_i$ $(e_{n+1},\cdots,e_m)=(-d_{n+1},\cdots,-d_m)$ $||e||_2^2=d_{n+1}^2+\cdots+d_m^2$ 。

通过QR分解实现最小二乘

给定m×n不一致系统
$A x = b$
找出完全QR分解A=QR，令
$\hat R=R$ 的上n×n子矩阵
$\hat d=d=Q^Tb$ 的上面的n个元素
$\hat R\bar x=\hat d$ $\bar x$

相比于法线方程求解最小二乘，这个更精确。，可以更好处理病态问题。

改进的格拉姆-施密特正交

改进的格拉姆-施密特正交

$A_j$ （j=1，...，n）为线性无关向量。
for j = 1，2，...，n
$y=A_j$
for i = 1，2，...，j-1
$r_{ij}=q_i^Ty$
$y=y-r_{ij}q_i$
end
$r_{jj}=||y||_2$
$q_j=y/r_{jj}$
end

$A_j$ $A_j$ $q_2$ $A_j$ $q_1$ $A_j$ $q_2$ 上。

豪斯霍尔德反射子

尽管改进的格拉姆-施密特正交是计算矩阵的QR分解的有效方式，但是这不是最好方式。另一种方法是使用豪斯霍尔德反射，这种方式需要更少的计算，同时在舍入误差放大的意义上来讲这种方法也更稳定。

豪斯霍尔德反射子是正交矩阵，通过m-1维平面反射m维向量。这意味着当每个向量乘上矩阵后，长度保持不变，使得豪斯霍尔德反射被称为移动向量的完美形式。给定一个向量x，我们要重新找出一个相同长度的向量w，计算豪斯霍尔德反射得出矩阵H满足Hx=w。

如图所示，原始方法非常清晰。画出m-1维平面二分x和w，并和连接它们的向量垂直，然后通过该平面反射所有向量。

引理 $||x||_2=||w||_2$ ，则w-x和w+x正交。

定义向量v=w-x，考虑投影矩阵

P = \frac{v v^{T}}{v^{T} v}

投影矩阵 $P^2=P$ 的矩阵。并且P满足Pv=v。几何上说，对于任何向量u，Pu是u在v上的投影。从上图可以看到，如果从x中两次减去Px，就可以得到w。令H=I-2P，则

\begin{aligned} H x = x - 2 P x = w - v - \frac{2 v v^{T} x}{v^{T} v} = w - v - \frac{v v^{T} x}{v^{T} v} - \frac{v v^{T} (w - v)}{v^{T} v} \\ = w - \frac{v v^{T} v}{v^{T} v} - \frac{v v^{T} x}{v^{T} v} - \frac{v v^{T} (w - v)}{v^{T} v} = w - \frac{v v^{T} w}{v^{T} v} - \frac{v v^{T} (w - v)}{v^{T} v} \\ = w - \frac{v v^{T} (2 w - v)}{v^{T} v} = w - \frac{v {(w - x)}^{T} (w + x)}{v^{T} v} = w \end{aligned}

矩阵H被称为豪斯霍尔德反射子，注意到H是对称并正交的矩阵。

定理（豪斯霍尔德反射子） $||x||_2=||w||_2$ $H=I-2vv^T/v^Tv$ 是对称正交矩阵，并且Hx=w。

使用豪斯霍尔德反射子推导出QR分解的一个新方式，使用反射子的目的是将列向量x移动到坐标轴，并以此将0放在矩阵中。

$x_1$ $\pm$ ${\left\| {{x_1}} \right\|_2}$ $H_1x=w$ $H_1$ 得到

H_{1} A = H_{1} [\begin{matrix} \times & \times & \times \\ \times & \times & \times \\ \times & \times & \times \\ \times & \times & \times \end{matrix}] = [\begin{matrix} \times & \times & \times \\ 0 & \times & \times \\ 0 & \times & \times \\ 0 & \times & \times \end{matrix}]

我们已经在A中引入一些0，我们希望继续这种方式直到A变为上三角；然后我们将得到QR分解中的R。对于一个3列的矩阵，还需要做两次这样的操作，即

\begin{aligned} H_{3} H_{2} H_{1} A = R \\ A = {H_{1}}^{- 1} {H_{2}}^{- 1} {H_{3}}^{- 1} R = H_{1} H_{2} H_{3} R = Q R \end{aligned}

举例：使用豪斯霍尔德反射子找到矩阵A的QR分解

A = [\begin{matrix} 1 & - 4 \\ 2 & 3 \\ 2 & 2 \end{matrix}]

$H_1$ $w=[||x||_2,0,0]$ =[3,0,0]，v=w-x=[2,-2,-2]

H_{1} = I - 2 \frac{v v^{T}}{v^{T} v} = [\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix}] - \frac{2}{12} [\begin{matrix} 4 & - 4 & - 4 \\ - 4 & 4 & 4 \\ - 4 & 4 & 4 \end{matrix}] = [\begin{matrix} \frac{1}{3} & \frac{2}{3} & \frac{2}{3} \\ \frac{2}{3} & \frac{1}{3} & - \frac{2}{3} \\ \frac{2}{3} & - \frac{2}{3} & \frac{1}{3} \end{matrix}]

H_{1} A = [\begin{matrix} \frac{1}{3} & \frac{2}{3} & \frac{2}{3} \\ \frac{2}{3} & \frac{1}{3} & - \frac{2}{3} \\ \frac{2}{3} & - \frac{2}{3} & \frac{1}{3} \end{matrix}] [\begin{matrix} 1 & - 4 \\ 2 & 3 \\ 2 & 2 \end{matrix}] = [\begin{matrix} 3 & 2 \\ 0 & - 3 \\ 0 & - 4 \end{matrix}]

$\hat H_2$ $\hat x$ $\hat w$ =[5,0]，

\begin{aligned} {\hat{H}}_{2} = I - 2 \frac{v v^{T}}{v^{T} v} = [\begin{matrix} 1 & 0 \\ 0 & 1 \end{matrix}] - \frac{2}{80} [\begin{matrix} 64 & - 32 \\ - 32 & 16 \end{matrix}] = [\begin{matrix} - 0.6 & - 0.8 \\ - 0.8 & 0.6 \end{matrix}] \\ {\hat{H}}_{2} {\hat{A}}_{2} = [\begin{matrix} - 0.6 & - 0.8 \\ - 0.8 & 0.6 \end{matrix}] [\begin{matrix} - 3 \\ - 4 \end{matrix}] = [\begin{matrix} 5 \\ 0 \end{matrix}] \end{aligned}

因此，有

H_{2} H_{1} A = [\begin{matrix} 1 & 0 & 0 \\ 0 & - 0.6 & - 0.8 \\ 0 & - 0.8 & 0.6 \end{matrix}] [\begin{matrix} \frac{1}{3} & \frac{2}{3} & \frac{2}{3} \\ \frac{2}{3} & \frac{1}{3} & - \frac{2}{3} \\ \frac{2}{3} & - \frac{2}{3} & \frac{1}{3} \end{matrix}] [\begin{matrix} 1 & - 4 \\ 2 & 3 \\ 2 & 2 \end{matrix}] = [\begin{matrix} 3 & 2 \\ 0 & 5 \\ 0 & 0 \end{matrix}] = R

A = H_{1} H_{2} R = Q R

与格拉姆-施密特正交相比，QR分解的计算代价更低，可以更好地得到单位正交向量以及具有更低的内存需求。

广义最小余项（GMRES）方法

Krylov方法

$A^kr$ $r=b-Ax_0$ $A^kr$ 对于大的k倾向于一个共有方向，Krylov空间的基必须认真计算。找出Krylov空间精确的基需要正交化计算方法，例如格拉姆-施密特或者豪斯霍尔德反射子方法。

$x_0$ $A^kr$ $x_0$ 中。

广义最小余项方法

$x_0$ =初始估计
$r=b-Ax_0$
$q_1=r/||r||_2$
for k = 1,2,...,m
$y=Aq_k$
for j = 1,2,...,k
$h_{jk}=q_j^Ty$
$y=y-h_{jk}q_j$
end
$h_{k+1,k}=||y||_2$ $h_{k+1,k}=0$ ，跳过下一行，并在底端终止）
$q_{k+1}=y/h_{k+1,k}$
$||Hc_k-[||r||_2,0,0,...,0]^T||_2$ $c_k$
$x_k=Q_kc_k+x_0$
end

$x_k$ $h_{k+1,k}=0$ $x_k$ $c_k$ $Q_k$ $Q_k$ $Q_k$ $x_k$ $x_0$ ，这种方法被称为重启GMRES。

预条件GMRES

$M^{-1}Ax=M^{-1}b$ ，其中M是方程组中讨论的预条件子。

预条件GMRES

$x_0$ =初始估计
$r=M^{-1}(b-Ax_0)$
$q_1=r/||r||_2$
for k = 1,2,...,m
$w=M^{-1}Aq_k$
for j = 1,2,...,k
$h_{jk}=w^Tq_j$
$w=w-h_{jk}q_j$
end
$h_{k+1,k}=||w||_2$
$q_{k+1}=w/h_{k+1,k}$
$||Hc_k-[||r||_2,0,0,...,0]^T||_2$ $c_k$
$x_k=Q_kc_k+x_0$
end

非线性最小二乘

高斯-牛顿方法

$F(x)^T=(r^TDr)^T=(Dr)^Tr$ ，

D F (x)^{T} = D ((D r)^{T} r) = (D r)^{T} \cdot D r + \sum_{i = 1}^{m} r_{i} D c_{i}

$c_i$ $Dc_i=Hr_i$ $r_i$ 的2阶偏导数矩阵或称为海森矩阵如下：

H_{r_{i}} = [\begin{matrix} \frac{\partial^{2} r_{i}}{\partial x_{1} \partial x_{1}} & \dots & \frac{\partial^{2} r_{i}}{\partial x_{1} \partial x_{n}} \\ ⋮ & ⋮ \\ \frac{\partial^{2} r_{i}}{\partial x_{n} \partial x_{1}} & \dots & \frac{\partial^{2} r_{i}}{\partial x_{n} \partial x_{n}} \end{matrix}]

高斯牛顿方法

为了最小化
$r_{1} (x)^{2} + \dots + r_{m} (x)^{2}$
$x^0$ =初始向量，
for k = 0,1,2,...
$A=Dr(x^k)$
$A^TAv^k=-A^Tr(x^k)$
$x^{k+1}=x^k+v^k$
end

D r (x, y) = [\begin{matrix} \frac{\partial r}{\partial x_{1}} & \frac{\partial r}{\partial y_{1}} \\ ⋮ & ⋮ \\ \frac{\partial r}{\partial x_{n}} & \frac{\partial r}{\partial y_{n}} \end{matrix}]

具有非线性参数的模型

$D_r$ $r_i$ $c_j$ 偏导数矩阵

(D r)_{i j} = \frac{\partial r_{i}}{\partial c_{j}} = f_{c_{j}} (t_{i})

举例 $c_1e^{c_2t}$

定义余项

r = [\begin{matrix} c_{1} e^{c_{2} t_{1}} - y_{1} \\ \dots \\ c_{1} e^{c_{2} t_{m}} - y_{m} \end{matrix}]

$c_1$ $c_2$ 计算偏导数得到

\begin{aligned} r = [\begin{matrix} c_{1} e^{c_{2} t_{1}} - y_{1} \\ \dots \\ c_{1} e^{c_{2} t_{m}} - y_{m} \end{matrix}] \\ D r = [\begin{matrix} e^{c_{2} t_{1}} & c_{1} t_{1} e^{c_{2} t_{1}} \\ ⋮ & ⋮ \\ e^{c_{2} t_{m}} & c_{1} t_{m} e^{c_{2} t_{m}} \end{matrix}] \end{aligned}

最小二乘问题中的非线性带来额外的挑战。法线方程以及QR方法只要系数矩阵A满秩都可以找到唯一解。而对于非线性问题的高斯-牛顿迭代可能收敛到多个极小平方误差中的一个。尽可能使用初始向量的合理近似，有助于收敛到绝对极小。

Levenberg-Marquardt方法

非线性最小二乘法遇到病态系数矩阵时，会面临巨大的挑战，Levenberg-Marquardt方法使用“正则化项”部分修复这个问题。

Levenberg-Marquardt方法

最小化
$r_{1} (x)^{2} + \dots + r_{m} (x)^{2}$
$x^0$ $\lambda$ =常数
for k = 0,1,2,...
$A=Dr(x^k)$
$(A^TA+\lambda diag(A^TA))v^k=-A^Tr(x^k)$
$x^{k+1}=x^k+v^k$
end

$A^TA$ $x_0$ 开始，并实现收敛。

数值微分和积分

计算微积分解决的主要问题是计算函数的导数和积分。

数值微分

有限差分公式

二点前向差分公式

f^{'} (x) = \frac{f (x + h) - f (x)}{h} - \frac{h}{2} f^{″} (c)

$- {h \over 2}f''(c)$ 是误差项。

定理（推广中值定理） $x_1$ $x_n$ $a_1$ $a_n$ >0，则在a和b之间存在数字c满足

(a_{1} + \dots + a_{n}) f (c) = a_{1} f (x_{1}) + \dots + a_{n} f (x_{n})

三点中心差分公式

f^{'} (x) = \frac{f (x + h) - f (x - h)}{2 h} - \frac{h^{2}}{6} f^{‴} (c)

其中x-h<c<x+h。

二阶导数的三点中心差分公式

f^{″} (x) = \frac{f (x - h) - 2 f (x) + f (x + h)}{h^{2}} - \frac{h^{2}}{12} f^{(4)} (c)

其中x-h<c<x+h。

当h->0时，误差会先减小再增大，后增大的原因是计算时有效数字的丢失。对于三点中心差分公式而言，f'(x)的及其近似误差的绝对值上界是

E (h) \equiv \frac{h^{2}}{6} f^{‴} (c) + \frac{ε_{m a c h}}{h}

外推

假设我们有n阶公式F(h)近似一个给定量Q，这个阶数意味着

Q \approx F (h) + K h^{″}

其中K大约是我们感兴趣的h区间上的一个常数，相关的例子是

f^{'} (x) = \frac{f (x + h) - f (x - h)}{2 h} - \frac{f^{‴} (c_{h})}{6} h^{2}

$c_h$ $f'''(c_h)$ $f'''(x)$ 相差不会太大。这样我们就可以把n阶公式变换成更高阶的公式，我们期望

Q - F (h / 2) \approx \frac{1}{2^{n}} (Q - F (h))

n阶公式的外推

Q \approx \frac{2^{n} F (h / 2) - F (h)}{2^{n} - 1}

举例：应用外推公式

\begin{aligned} F_{4} (x) = \frac{2^{2} F_{2} (h / 2) - F_{2} (h)}{2^{2} - 1} \\ = [4 \frac{f (x + h / 2) - f (x - h / 2)}{h} - \frac{f (x + h) - f (x - h)}{2 h}] / 3 \\ = \frac{f (x - h) - 8 f (x - h / 2) + 8 f (x + h / 2) - f (x + h)}{6 h} \end{aligned}

数值积分的牛顿-科特斯公式

梯形法则

\int_{x_{0}}^{x_{1}} f (x) d x = \frac{h}{2} (y_{0} + y_{1}) - \frac{h^{3}}{12} f^{″} (c)

$h=x_1-x_0$ $x_0$ $x_1$ 之间。

辛普森法则

\int_{x_{0}}^{x_{2}} f (x) d x = \frac{h}{3} (y_{0} + 4 y_{1} + y_{2}) - \frac{h^{5}}{90} f^{(4)} (c)

$h=x_2-x_1=x_1-x_0$ $x_0$ $x_2$ 之间。

定义：数值积分方法的精度是最大的整数k，使用该积分方法可以得到所有k阶或者更低阶多项式积分的精确值。

梯形法则的精度是1，辛普森法则的精度是3。

复合牛顿-科特斯公式

由于积分在区间的所有子区间上具有可加性，我们可以通过除法把整个区间变为很多小区间再计算积分，在每个小区间上使用法则，然后再求和。这种策略被称为复合数值积分。

复合梯形法则

\int_{a}^{b} f (x) d x = \frac{h}{2} (y_{0} + y_{m} + 2 \sum_{i = 1}^{m - 1} y_{i}) - \frac{(b - a) h^{2}}{12} f^{″} (c)

其中h=(b-a)/m，c在a和b之间。

复合辛普森公式

\int_{a}^{b} f (x) d x = \frac{h}{3} (y_{0} + y_{2 m} + 4 \sum_{i = 1}^{m} y_{2 i - 1} + 2 \sum_{i = 1}^{m - 1} y_{2 i}) - \frac{(b - a) h^{4}}{180} f^{(4)} (c)

其中c在a和b之间。

开牛顿-科特斯方法

中点法则

\int_{x_{0}}^{x_{1}} f (x) d x = h f (w) + \frac{h^{3}}{24} f^{″} (c)

$h=(x_1-x_0)$ $x_0+h/2$ $x_0$ $x_1$ 之间。

中点法则也可用于减少所需函数求值的次数与梯形法则相比，闭牛顿-科特斯方法具有相同的阶数，只需要一次函数求值，而不是两次。而且误差项是梯形法则误差项大小的一半。

复合中点法则

\int_{x_{0}}^{x_{1}} f (x) d x = h \sum_{i = 1}^{m} f (w_{i}) + \frac{(b - a) h^{2}}{24} f^{″} (c)

$w_i$ 是[a,b]中m个相等子区间的中点。

龙贝格积分

龙贝格积分是对复合梯形法则应用外推的结果，给定近似M的法则N(h)，该法则依赖步长h，如果知道法则的阶，则可以对法则进行外推。

龙贝格积分

$R_{11}=(b-a){f(a)+f(b)\over 2}$
for j = 2,3,...n
$h_j={b-a\over 2^{j-1}}$
${R_{j1}} = {1 \over 2}{R_{j - 1,1}} + {h_j}\sum\limits_{i = 1}^{{2^{j - 2}}} {f(a + (2i - 1){h_j})}$
for k = 2,...,j
${R_{jk}} = {{{4^{k - 1}}{R_{j,k - 1}} - {R_{j - 1,k - 1}}} \over {{4^{k - 1}} - 1}}$
end
end

$R_{jj}$ $R_{jj}$ 差异小于当前的容差。

自适应积分

到目前我们学到的积分近似方法都使用相等的步长。一般来说，小步长可以提高精度变化剧烈的函数将需要更多步，因而带来更多的计算时间，这是由于需要更小的步长跟踪函数所有的变化。函数可能在定义域某些部分变化剧烈，而在另外一些部分变化缓慢，满足前一部分的容差的步长在后一部分就显得过于精细了。

通过使用积分误差公式，可以在运算中推出一个标准，其步长对于特定的子空间适合。这种方法称为自适应积分。

自适应积分

$\int_a^bf(x)dx$
c = (a+b)/2
$S_{[a,b]}$ =(b-a)(f(a)+f(b))/2
$S_{[a,b]}-S_{[a,c]}-S_{[c,b]}<3\cdot TOL\cdot ({b-a \over b_{orig}-a_{orig}})$
$S_{[a,c]}+S_{[c,b]}$ 作为区间[a,b]上的近似
else
对于区间[a,c]和[c,b]递归重复上面的步骤
end


xxxxxxxxxx
% 程序 5 . 2 自适应积分
% 计算定积分的近似
% 输入： MATLAB函数f，区间 [ ao, bo]
% 容差 tol0
% 输出：近似定积分
function int=adapquad(f, a0, b0, tol0)
int = 0; n=1; a(1) = a0; b(1)=b0; tol(1)=tol0; 
app(1)=trap(f,a,b);                         % trap求梯形积分
while n >0                                  % n是当前列表结束的位置。
    c = (a(n) + b(n)) / 2; oldapp = app(n);
    app(n) = trap(f, a(n), c);
    app(n+1) = trap(f, c, b(n));
    if abs(oldapp - (app(n) + app(n+1))) < 3 * tol(n)
        int = int + app(n) + app(n+1);      % 成功
        n = n-1;                            % 该区间操作完成
    else                                    % 分成两个子区间
        b(n+1) = b(n); b(n) = c;            % 设置新区间
        a(n+1) = c;
        tol(n) = tol(n) / 2; tol(n+1) = tol(n);
        n = n + 1;                          % 到列表尾端，重复
    end
end

其中tol0是自己设置的参数，如设置为0.005。

高斯积分

定义 $p_0$ $p_n$ }在区间[a,b]正交，当

\int_{a}^{b} p_{j} (x) p_{k} (x) d x = {\begin{matrix} 0 j \neq k \\ \neq 0 j = k \end{matrix}

定理 $p_0$ $p_n$ $\deg p_i=i$ $p_0$ $p_n$ }是在区间[a,b]由区间[a,b]由至多n个多项式所张的向量空间的基。

定理 $p_0$ $p_n$ $\deg p_i=i$ $p_i$ 在区间(a,b)上有i个不同的根。

高斯积分

\int_{- 1}^{1} f (x) d x \approx \sum_{i = 1}^{n} c_{i} f (x_{i})

其中

c_{i} = \int_{- 1}^{1} L_{i} (x) d x, i = 1, \dots, n

L_{i} = \frac{(x - x_{1}) \dots \overset{―}{(x - x_{i})} \dots (x - x_{n})}{(x_{i} - x_{1}) \dots \overset{―}{(x_{i} - x_{i})} \dots (x_{i} - x_{n})}

$c_i$ 可以查表获得(包括n=1,2,3,4)

定理：高斯积分方法，在区间[—1,1]上使用n阶勒让德多项式，具有2n-1阶精度。

一般区间[a,b]的积分可以通过替代t=(2x-a-b)/(b-a)得到

\int_{a}^{b} f (x) d x = \int_{- 1}^{1} f (\frac{(b - a) t + b + a}{2}) \frac{b - a}{2} d t

常微分方程

初值问题

欧拉方法

$t_0$ $y_0$ $t_1$ $y_1$ )重新计算斜率，根据新的斜率继续移动，重复这个过程。

欧拉方法

$w_0=y_0$
$w_{i+1}=w_i+hf(t_i,w_i)$

h是步长，如取0.2，y'=f(t,w)。

举例：用欧拉方法求解如下问题

{\begin{matrix} y^{'} = t y + t^{3} \\ y (0) = y_{0} \\ t \in [0, 1] \end{matrix}

欧拉方法则对应如下迭代，h设置为0.1

\begin{aligned} w_{0} = 1 \\ w_{i + 1} = w_{i} + h (t_{i} w_{i} + t_{i}^{3}) \end{aligned}

解的存在性、唯一性和连续性

定义利普希茨常数 $y_1$ $y_2$ )满足

| f (t, y_{1}) - f (t, y_{2}) | \leq L | y_{1} - y_{2} |

函数f(t,y)相对于变量y在S上利普西茨连续。

定理 $y_a$ <β，函数对于变量y是利普希茨连续，则在a与 b之间存在c，使得初值问题

{\begin{matrix} y^{'} = f (t, y) \\ y (a) = y_{a} \\ t \in [a, c] \end{matrix}

有唯一解y(t)。而且，如果f在[a,b]×(-∞，∞)是利普西茨连续，则它在[a,b]上存在唯一解。

定理：假设f(t,y)在集合S=[a,b]×[α,β]上关于y是连续的，如果Y(t)和Z(t)是微分方程

y^{'} = f (t, y)

在S上的解，分别具有初值条件Y(z)和Z(a)，则

| Y (t) - Z (t) | \leq e^{L (t - a)} | Y (a) - Z (a) |

一阶线性方程

一组容易求解的特定的常微分方程提供了ODE求解问题中具有指导意义的例子，这组方程是一阶方程，其右侧是关于y的线性函数。考虑初值问题

{\begin{matrix} y^{'} = g (t) y + h (t) \\ y (a) = y_{a} \\ t \in [a, b] \end{matrix}

$L=\max_{[a,b]}g(t)$ 作为利普希茨常数可以使用技巧找到解，即对方程乘上“积分因子”。

${e^{ - \int {g(t)} dt}}$ ，在方程两侧同时乘这个因子得到

\begin{aligned} (y^{'} - g (t) y) e^{- \int g (t) d t} = e^{- \int g (t) d t} h (t) \\ (y e^{- \int g (t) d t})^{'} = e^{- \int g (t) d t} h (t) \\ y e^{- \int g (t) d t} = \int e^{- \int g (t) d t} h (t) d t \\ y (t) = e^{\int g (t) d t} \int e^{- \int g (t) d t} h (t) d t \end{aligned}

IVP求解器的分析

显式梯形方法

显式梯形方法

$w_0=y_0$
${w_{i + 1}} = {w_i} + {h \over 2}\left( {f({t_i},{w_i}) + f({t_i} + h,{w_i} + hf({t_i},{w_i}))} \right)$

泰勒方法

泰勒方法基本思想是直接利用泰勒展开，假设解y(t)是(k+1)阶连续可微的函数。给定在解曲线上的当前点(t,y(t))，目标是对于某个步长h，使用微分方程的信息，用y(t)来表达y(t+h)。

K阶泰勒方法

$w_0=y_0$
${w_{i + 1}} = {w_i} + hf({t_i},{w_i}) + {{{h^2}} \over 2}f'({t_i},{w_i}) + \cdots + {{{h^k}} \over {k!}}{f^{(k - 1)}}({t_i},{w_i})$

举例：对于一阶线性方程确定二阶泰勒方法

{\begin{matrix} y^{'} = t y + t^{3} \\ y (0) = y_{0} \end{matrix}

$f(t,y)=ty+t^3$ ，

f^{'} (t, y) = f_{t} + f_{y} f = y + 3 t^{2} + t (t y + t^{3})

所以二阶泰勒方法的迭代公式为

w_{i + 1} = w_{i} + h (t_{i} w_{i} + t_{i}^{3}) + \frac{h^{2}}{2} (w_{i} + 3 {t_{i}}^{2} + t_{i} (t_{i} w_{i} + {t_{i}}^{3}))

常微分方程组

高阶方程

单个的高阶微分方程可以转化为一个方程组，令

y^{(n)} = f (t, y, y^{'}, y^{″}, \dots, y^{(n - 1)})

是n阶常微分方程，定义新的变量

\begin{aligned} y_{1} = y \\ y_{2} = y^{'} \\ y_{3} = y^{″} \\ ⋮ \\ y_{n} = y^{(n - 1)} \end{aligned}

注意到原始的常微分方程写成

y_{n}^{'} = f (t, y_{1}, y_{2}, \dots, y_{n})

两者放在一起

\begin{aligned} {y_{1}}^{'} = y_{2} \\ {y_{2}}^{'} = y_{3} \\ {y_{3}}^{'} = y_{4} \\ ⋮ \\ {y_{n - 1}}^{'} = y_{n} \\ y_{n}^{'} = f (t, y_{1}, y_{2}, \dots, y_{n}) \end{aligned}

这些方程将n阶的微分方程转化为一阶微分方程组，该方程组可以使用欧拉或者梯形方法求解。

龙格-库塔方法和应用

龙格-库塔家族

我们已经知道一阶的欧拉方法以及二阶的梯形方法。除了梯形方法外，还有其他龙格-库塔类型的二阶方法。

中点方法

$w_0=y_0$
${w_{i + 1}} = {w_i} + hf\left( {{t_i} + {h \over 2},{w_i} + {h \over 2}f({t_i},{w_i})} \right)$

4阶龙格-库塔方法（RK4）

${w_{i + 1}} = {w_i} + {h \over 6}\left( {{s_1} + 2{s_2} + 2{s_3} + {s_4}} \right)$
其中
$\begin{aligned} s_{1} = f (t_{i}, w_{i}) \\ s_{2} = f (t_{i} + \frac{h}{2}, w_{i} + \frac{h}{2} s_{1}) \\ s_{3} = f (t_{i} + \frac{h}{2}, w_{i} + \frac{h}{2} s_{2}) \\ s_{4} = f (t_{i} + h, w_{i} + h s_{3}) \end{aligned}$

可变步长方法

直到现在，步长h在ODE求解器实现中一直是一个常数。但是，并没有原因说h在求解过程中不能改变。我们希望改变步长的一个原因是，问题的解会在缓慢变化的周期和迅速变化的周期之间移动。将固定步长变得足够小使其精确跟踪快速变化，意味着解的其他部分求解都会是令人难以忍受的缓慢。

龙格-库塔嵌入对

可变步长方法的关键思想是检测当前步生成的误差。用户设置的容差在当前步必须能够满足．然后设计方法(1)如果超出容差，必须拒绝误差并将步长减小，或者(2)如果满足容差，接受步长并选择对于下一步适合的步长h。关键是近似在每步中产生的误差。首先假设我们已经找到这样的方式并解释如何改变步长。

有RK2/3（龙格-库塔2阶/3阶嵌入对）、RK4/5，不过多赘述。

隐式方法和刚性方程

后向欧拉方法

$w_0=y_0$
$w_{i+1}=w_i+hf(t_{i+1},w_{i+1})$

举例：对初值问题使用后向欧拉方法

{\begin{matrix} y^{'} = y + 8 y^{2} - 9 y^{3} \\ y (0) = 1 / 2 \\ t \in [0, 3] \end{matrix}

使用后向欧拉方法

w_{i + 1} = w_{i} + h f (t_{i + 1}, w_{i + 1}) = w_{i} + h (w_{i + 1} + 8 w_{i + 1}^{2} - 9 w_{i + 1}^{3})

$z=w_{i+1}$ $z=w_i+h(z+8z^2-9z^3)$ ，或者对于未知的z求解

9 h z^{3} - 8 h z^{2} + (1 - h) z - w_{i} = 0

可以使用牛顿方法进行估计

z_{n e w} = z - \frac{9 h z^{3} - 8 h z^{2} + (1 - h) z - w_{i}}{27 h z^{2} - 16 h z + 1 - h}

$z_{new}$ $z_{new}-z$ 小于容差（比近似微分方程解的误差还要小）。

多步方法

$w_{i+1}$ $w_i$ $w_i$ 之外的知识帮助下一步的求解。这将带来OD E求解器的阶数和单步方法的阶数相同，但是大量必要的计算会被替换为求解过程中已经计算的值的插值。

构造多步方法

Adams-Bashforth两步方法

w_{i + 1} = w_{i} + h [\frac{3}{2} f (t_{i}, w_{i}) - \frac{1}{2} f (t_{i - 1}, w_{i - 1})]

显式多步方法

Adams-Bashforth三步方法（三阶）

w_{i + 1} = w_{i} + \frac{h}{12} [23 f_{i} - 16 f_{i - 1} + 5 f_{i - 2}]

Adams-Bashforth四步方法（四阶）

w_{i + 1} = w_{i} + \frac{h}{24} [55 f_{i} - 59 f_{i - 1} + 37 f_{i - 2} - 9 f_{i - 3}]

隐式多步方法

隐式梯形方法（二阶）

w_{i + 1} = w_{i} + \frac{h}{2} [f_{i + 1} + f_{i}]

Adams-Moulton两步方法（三阶）

w_{i + 1} = w_{i} + \frac{h}{12} [5 f_{i + 1} + 8 f_{i} - f_{i - 1}]

Milne-Simpson方法

w_{i + 1} = w_{i} + \frac{h}{3} [f_{i + 1} + 4 f_{i} + f_{i - 1}]

Adams-Moulton三步方法（四阶）

w_{i + 1} = w_{i} + \frac{h}{24} [9 f_{i + 1} + 19 f_{i} - 5 f_{i - 1} + f_{i - 2}]

Adams-Moulton四步方法（五阶）

w_{i + 1} = w_{i} + \frac{h}{720} [251 f_{i + 1} + 646 f_{i} - 246 f_{i - 1} + 106 f_{i - 2} - 19 f_{i - 3}]

边值问题

有限差分方法

线性边值问题

令y(t)是一个至少四阶连续的函数，可以得到一阶导数的离散近似为

y^{'} (t) = \frac{y (t + h) - y (t - h)}{2 h} - \frac{h^{2}}{6} y^{‴} (c)

二阶导数的近似为

\begin{aligned} y^{'} (t) = \frac{y (t + h) - y (t - h)}{2 h} - \frac{h^{2}}{6} y^{‴} (c) \\ y^{″} (t) = \frac{y (t + h) - 2 y (t) + y (t - h)}{h^{2}} + \frac{h^{2}}{12} y^{⁗} (c) \end{aligned}

$y_i$ $w_i$ ，边界条件则在需要的时候代入方程组中。

举例：使用有限差分求解

{\begin{matrix} y^{″} = 4 y \\ y (0) = 1 \\ y (1) = 3 \end{matrix}

$t_i$ 处的有限差分形式如下

\frac{w_{i + 1} - 2 w_{i} + w_{i - 1}}{h^{2}} - 4 w_{i} = 0

或者等价地

w_{i - 1} + (- 4 h^{2} - 2) w_{i} + w_{i + 1} = 0

$w_0=1$ $w_4=3$ ，给出方程

{\begin{matrix} 1 + (- 4 h^{2} - 2) w_{1} + w_{2} = 0 \\ w_{1} + (- 4 h^{2} - 2) w_{2} + w_{3} = 0 \\ w_{2} + (- 4 h^{2} - 2) w_{3} + 3 = 0 \end{matrix}

代入h得到三对角线矩阵方程

[\begin{matrix} - \frac{9}{4} & 1 & 0 \\ 1 & - \frac{9}{4} & 1 \\ 0 & 1 & - \frac{9}{4} \end{matrix}] [\begin{matrix} w_{1} \\ w_{2} \\ w_{3} \end{matrix}] = [\begin{matrix} - 1 \\ 0 \\ - 3 \end{matrix}]

非线性边值问题

对于非线性边值问题，使用多变量牛顿方法进行迭代，迭代公式为

w^{k + 1} = w^{k} - D F (w^{k})^{- 1} F (w^{k})

$DF({w^k})\Delta w = - F({w^k})$ $\Delta w = {w^{k + 1}} - {w^k}$ 。

举例：求解非线性边值问题

{\begin{matrix} y^{″} = y - y^{2} \\ y (0) = 1 \\ y (1) = 4 \end{matrix}

离散化后方程为

w_{i - 1} - (2 + h^{2}) w_{i} + h^{2} w_{i}^{2} + w_{i + 1} = 0

$w_0=y_a=1$ $w_{n+1}=y_b=4$ 。

函数F(w)如下

F [\begin{matrix} w_{1} \\ w_{2} \\ ⋮ \\ w_{n - 1} \\ w_{n} \end{matrix}] = [\begin{matrix} \begin{matrix} y_{a} - (2 + h^{2}) w_{1} + h^{2} w_{1}^{2} + w_{2} \\ w_{1} - (2 + h^{2}) w_{2} + h^{2} w_{2}^{2} + w_{3} \end{matrix} \\ ⋮ \\ w_{n - 2} - (2 + h^{2}) w_{n - 1} + h^{2} w_{n - 1}^{2} + w_{n} \\ w_{n - 1} - (2 + h^{2}) w_{n} + h^{2} w_{n}^{2} + y_{b} \end{matrix}]

排列与有限元方法

和有限差分方法相似，排列和有限元方法背后的思想是将边值问题消减为一组可求解的代数方程。但是，并不是通过使用有限差分替换微分方程中的导数进行离散化，而是对于解给出函数形式，其对应的参数使用该方法拟合。

${\phi _1}(t)$ ${\phi _n}(t)$ ，它们可能是多项式、三角函数、样条，或者其他简单函数，然后考虑可能的解

y (t) = c_{1} ϕ_{1} (t) + \dots + c_{n} ϕ_{n} (t)

$c_i$ 的值。

排列

考虑BVP问题

{\begin{matrix} y^{″} = f (t, y, y^{'}) \\ y (a) = y_{a} \\ y (b) = y_{b} \end{matrix}

选择n个点，从边值点a开始，在b点结束，即

a = t_{1} < t_{2} < \dots < t_{n} = b

$c_1$ $c_n$ $t_i$ ${t_i} = a + {{i - 1} \over {n - 1}}\left( {b - a} \right)$ 。

$\phi_j(t)=t^{j-1}$ ，其中1≤j≤n。解的形式如下

y (t) = \sum_{j = 1}^{n} c_{j} ϕ_{j} (t) = \sum_{j = 1}^{n} c_{j} t^{j - 1}

$c_1$ $c_n$ 的n个方程，其中第一个和最后一个方程是边值条件

\begin{aligned} \sum_{j = 1}^{n} c_{j} a^{j - 1} = y (a) \\ \sum_{j = 1}^{n} c_{j} b^{j - 1} = y (b) \end{aligned}

$t_i$ $y(t)=\sum_{j=1}^nc_jt^{j-1}$ 得到

\sum_{j = 1}^{n} (j - 1) (j - 2) c_{j} t^{j - 3} = f (t, \sum_{j = 1}^{n} c_{j} t^{j - 1}, \sum_{j = 1}^{n} (j - 1) c_{j} t^{j - 2})

$t_i$ $c_i$ $c_i$ 的方程组也是线性的，并可以求解。

在排列方法中选择三角函数作为基带来了傅里叶分析以及谱方法，这在边值方程和偏微分方程中都大量应用。这是一个“全局的方法，其中对于t的一个非常大的范围中，基函数非0，但是具有很好的正交性质。

有限元以及Galerkin方法

$c_i$ 的一个不同的方程组。

考虑一个BVP问题

{\begin{matrix} y^{″} = f (t, y, y^{'}) \\ y (a) = y_{a} \\ y (b) = y_{b} \end{matrix}

解为y，余项r=y"-f，需要使得微分方程两侧的差异尽可能小。和最小二乘方法类似，这可以通过选择y使得余项与潜在解的向量空间正交得到。

对于区间[a,b]，定义平方可积分的函数对应的向量空间

L^{2} [a, b] = {区 间 [a, b] 上 的 函 数 y (t) | \int_{a}^{b} y {(t)}^{2} d t 存 在 并 有 穷}

$L^2$ 函数空间具有内积

⟨ y_{1}, y_{2} ⟩ = \int_{a}^{b} y_{1} (t) y_{2} (t) d t

$<y_1,y_2>=0$ $y_1$ $y_2$ $L^2[a,b]$ 正交 $L^2[a,b]$ $L^2[a,b]$ $L^2$ $\phi_0(t)$ $\phi_{n+1}(t)$ 。

$L^2$ $\int_a^b {\left( {y'' - f} \right){\phi _i}(t)} dt=0$ ，或

\int_{a}^{b} y^{″} ϕ_{i} (t) d t = \int_{a}^{b} f (t, y, y^{'}) ϕ_{i} (t) d t

其中0≤i≤n+1，该形式被称为边值问题的弱形式。

Galerkin方法的第二个思想是部分使用积分消去二阶导数，注意到

\begin{aligned} \int_{a}^{b} y^{″} ϕ_{i} (t) d t = {ϕ_{i} (t) y^{'} (t) |}_{a}^{b} - \int_{a}^{b} y^{'} (t) {ϕ_{i}}^{'} (t) d t \\ = ϕ_{i} (b) y^{'} (b) - ϕ_{i} (a) y^{'} (a) - \int_{a}^{b} y^{'} (t) {ϕ_{i}}^{'} (t) d t \end{aligned}

可以得到

\int_{a}^{b} f (t, y, y^{'}) ϕ_{i} (t) d t = ϕ_{i} (b) y^{'} (b) - ϕ_{i} (a) y^{'} (a) - \int_{a}^{b} y^{'} (t) {ϕ_{i}}^{'} (t) d t

$c_i$

y (t) = \sum_{i = 0}^{n + 1} c_{i} ϕ_{i} (t)

$\phi_i(t)$ ${t_0} < {t_1} < \cdots < {t_{n-1}} < {t_n}$ 开始，对于i=1，...，n定义

ϕ_{i} (t) = {\begin{matrix} \frac{t - t_{i - 1}}{t_{i} - t_{i - 1}} t_{i - 1} < t \leq t_{i} \\ \frac{t_{i + 1} - t}{t_{i + 1} - t_{i}} t_{i} < t \leq t_{i + 1} \\ 0 o t h e r w i s e \end{matrix}

同时定义

ϕ_{0} (t) = {\begin{matrix} \frac{t_{1} - t}{t_{1} - t_{0}} t_{0} \leq t < t_{1} \\ 0 o t h e r w i s e \end{matrix} ϕ_{n + 1} (t) = {\begin{matrix} \frac{t - t_{n}}{t_{n + 1} - t_{n}} t_{n} < t \leq t_{n + 1} \\ 0 o t h e r w i s e \end{matrix}

$\phi_i(t_j)$ $(t_i,c_i)$ ，定义分段线性B样条

S (t) = \sum_{i = 0}^{n + 1} c_{i} ϕ_{i} (t)

$S(t) = \sum\limits_{i = 0}^{n + 1} {{c_i}{\phi _i}(t_j)}=c_j$ 。

$c_i$ 通过排列方法得到

\begin{aligned} y (a) = \sum_{i = 0}^{n + 1} c_{i} ϕ_{i} (a) = c_{0} ϕ_{0} (a) = c_{0} \\ y (b) = \sum_{i = 0}^{n + 1} c_{i} ϕ_{i} (b) = c_{n + 1} ϕ_{n + 1} (b) = c_{n + 1} \end{aligned}

对于i=1，...，n，使用有限元方程

\int_{a}^{b} f (t, y, y^{'}) ϕ_{i} (t) d t + \int_{a}^{b} y^{'} (t) {ϕ_{i}}^{'} (t) d t = 0

$y(t)=\sum c_i\phi_i(t)$ $\phi$ 函数的性质可以求解。

偏微分方程

前向差分方法

定义偏微分方程（PDE）如下

{\begin{matrix} u_{t} = D u_{x x} a \leq x \leq b, t \geq 0 \\ u (x, 0) = f (x) a \leq x \leq b \\ u (a, t) = f (t) t \geq 0 \\ u (b, t) = r (t) t \geq 0 \end{matrix}

使用数值微分和积分的离散公式近似在x和t方向的导数。例如，使用中心差分公式，关于x的二阶导数如下

u_{x x} (x, t) \approx \frac{1}{h^{2}} (u (x + h, t) - 2 u (x, t) + u (x - h, t))

${h^2}{u_{xxxx}}({c_1},t)/12$ ，对于时间变量的一阶导数的前向差分公式得到

u_{t} (x, t) = \frac{1}{k} (u (x, t + k) - u (x, t))

$ku_{tt}(x,c_2)/2$ $x-h<c_1<x+h$ $t<c_2<t+h$ $(x_i,t_j)$ 处的偏微分方程中得到

\frac{D}{h^{2}} (w_{i + 1, j} - 2 w_{i j} + w_{i - 1, j}) \approx \frac{1}{k} (w_{i, j + 1} - w_{i j})

$\sigma={DK \over h^2}$ ）

\begin{aligned} w_{i, j + 1} = w_{i j} + \frac{D K}{h^{2}} (w_{i + 1, j} - 2 w_{i j} + w_{i - 1, j}) \\ = σ w_{i + 1, j} + (1 - 2 σ) w_{i j} + σ w_{i - 1, j} \end{aligned}

写成矩阵形式

[\begin{matrix} w_{1. j + 1} \\ ⋮ \\ w_{m, j + 1} \end{matrix}] = [\begin{matrix} 1 - 2 σ & σ & 0 & \dots & 0 \\ σ & 1 - 2 σ & σ & ⋱ & 0 \\ 0 & σ & 1 - 2 σ & ⋱ & 0 \\ ⋮ & ⋱ & ⋱ & ⋱ & σ \\ 0 & \dots & 0 & σ & 1 - 2 σ \end{matrix}] [\begin{matrix} w_{1. j} \\ ⋮ \\ w_{m, j} \end{matrix}] + σ [\begin{matrix} w_{0, j} \\ 0 \\ ⋮ \\ \begin{matrix} 0 \\ w_{m + 1, j} \end{matrix} \end{matrix}]

后向差分方法

$\sigma={DK \over h^2}$ ）

\begin{aligned} \frac{1}{k} (w_{i j} - w_{i, j - 1}) = \frac{D}{h^{2}} (w_{i + 1, j} - 2 w_{i j} + w_{i - 1, j}) \\ - σ w_{i + 1, j} + (1 + 2 σ) w_{i j} - σ w_{i - 1, j} = w_{i, j - 1} \end{aligned}

矩阵形式为

[\begin{matrix} 1 + 2 σ & - σ & 0 & \dots & 0 \\ - σ & 1 + 2 σ & - σ & ⋱ & 0 \\ 0 & - σ & 1 + 2 σ & ⋱ & 0 \\ ⋮ & ⋱ & ⋱ & ⋱ & - σ \\ 0 & \dots & 0 & - σ & 1 + 2 σ \end{matrix}] [\begin{matrix} w_{1. j} \\ ⋮ \\ w_{m, j} \end{matrix}] = [\begin{matrix} w_{1. j - 1} \\ ⋮ \\ w_{m, j - 1} \end{matrix}] + σ [\begin{matrix} w_{0, j} \\ 0 \\ ⋮ \\ \begin{matrix} 0 \\ w_{m + 1, j} \end{matrix} \end{matrix}]

Crank-Nicolson方法

Crank-Nicolson对于时间导数使用后向差分公式，并均匀组合前向差分近似和后向差分近似。

$u_t$

\frac{1}{k} (w_{i j} - w_{i, j - 1})

$u_{xx}$

\frac{1}{2} (\frac{w_{i + 1, j} - 2 w_{i j} + w_{i - 1, j}}{h^{2}}) + \frac{1}{2} (\frac{w_{i + 1, j - 1} - 2 w_{i, j - 1} + w_{i - 1, j - 1}}{h^{2}})

可以重新组织热方程近似

2 w_{i j} - 2 w_{i, j - 1} = σ [w_{i + 1, j} - 2 w_{i j} + w_{i - 1, j} + w_{i + 1, j - 1} - 2 w_{i, j - 1} + w_{i - 1, j - 1}]

或者

- σ w_{i - 1, j} + (2 + 2 σ) w_{i j} - σ w_{i + 1, j} = σ w_{i - 1, j - 1} + (2 - 2 σ) w_{i, j - 1} + σ w_{i + 1, j - 1}

本章中有各种偏微分方程（抛物线方程，双曲线方程，椭圆方程，非线性偏微分方程）的解法。

随机数和应用

随机数

伪随机数

定义 $u_i$ 为生成的[0,1]随机数）

\begin{matrix} x_{i} = a x_{i - 1} + b (mod m) \\ u_{i} = \frac{x_{i}}{m} \end{matrix}

其中a为乘子，b为偏移，m为模数。

最小标准随机数生成器

$x_i=ax_{i-1}(\bmod m)$
$u_i={x_i\over m}$
$m=2^{31}-1$ $a=7^5=16807$ $b=0$

按照现在计算机每秒指令数的增加，这个数目显得太少了。

randu生成器

$x_i=ax_{i-1}(\bmod m)$
$u_i={x_i\over m}$
$a=65539=2^{16}+3$ $m=2^{31}$

不满足随机数的独立性要求。

最新版的matlab已经不再使用LCG作为随机数生成器，从matlab5开始，使用的是延迟斐波那契生成器。

指数和正态随机数

$p(x)=ae^{-ax},a>0$ ，使用均匀分布随机数，可以很容易地生成指数随机数，累计分布函数为

P (x) = P r o b (V \leq x) = \int_{0}^{x} p (x) d x = 1 - e^{- a x}

算法的核心思想是要找到指数随机变量的值，使得概率Prob(V≤x)是[0,1]区间上的均匀分布，即给定一个满足均匀分布的变量u之后，要使得

u = P r o b (V \leq x) = 1 - e^{- a x}

求解x得到

x = \frac{- \ln (1 - u)}{a}

把输入的均匀随机数u转换为指数随机数。

$Q(x)=P^{-1}(x)$ 是对应的反函数。若U[0,1]是[0,1]区间上的均匀分布的随机数，则Q(U[0,1])就是满足分布P的随机数，剩下的问题就是如何高效地计算函数Q。

$p(x,y)=(1/2\pi)e^{-(x^2+y^2)/2}$ $p(r)=(1/2\pi)e^{-r^2/2}$ $r^2$ $r^2$

r^{2} = \frac{- \ln (1 - u_{1})}{1 / 2}

$u_1$ 是均匀分布随机数，进而生成两个满足标准正态分布且完全独立的随机数

\begin{aligned} n_{1} = r \cos 2 π u_{2} = \sqrt{- 2 \ln (1 - u_{1})} \cos 2 π u_{2} \\ n_{2} = r \sin 2 π u_{2} = \sqrt{- 2 \ln (1 - u_{1})} \sin 2 π u_{2} \end{aligned}

这种算法称为Box-Muller方法，此外还有一种更加高效的方法

n_{1} = x_{1} \sqrt{\frac{- 2 \ln (u_{1})}{u_{1}}} n_{2} = x_{2} \sqrt{\frac{- 2 \ln (u_{1})}{u_{1}}}

$x_1^2+x_2^2<1$ 就接受这个值，反之重新选择。

蒙特卡洛模拟

幂律和蒙特卡洛模拟

$E[(E-E(X))^2]$ ，而X标准差为方差的平方根。在测量过程中，误差的期望会随着n的增加而下降，并且符合如下的规律：

Ⅰ型或Ⅱ型蒙特卡罗模拟具有伪随机数字。

E r r o r \propto n^{- \frac{1}{2}}

拟随机数

拟随机数的理念是在条件允许的情况下，放弃对随机数独立性的要求。放弃独立性意味着拟随机数不再是随机的，而且也不是像伪随机数那样看上去是随机的。通过这种方式，可以使得蒙特卡罗方法能很快地收敛到正确解。拟随机数的设计是要使得生成的序列是自避(Self-avoiding)的，而不一定要保证独立性。自避定义为，在生成随机数序列的过程中，新的数会填充到比较稀疏的区域，而不会聚集到一起。

Halton序列 $b_kb_{k-1}\cdots b_2b_1$ $0.b_1b_2\cdots b_{k-1}b_k$ 。换句话说，整个生成过程就是先写出第i个整数，然后把位数反过来写，再放到小数点的另一侧，最后输出[0,1]区间上均匀分布的第i个随机数。

如p=3生成的前8个随机数

i	$(i)_3$	$(u_i)_3$	$u_i$	i	$(i)_3$	$(u_i)_3$	$u_i$
1	1	0.1	$0.\bar 3$	5	12	0.21	$0.\bar 7$
2	2	0.2	$0.\bar 6$	6	20	0.02	$0.\bar 2$
3	10	0.01	$0.\bar 1$	7	21	0.12	$0.\bar 5$
4	11	0.11	$0.\bar 4$	8	22	0.22	$0.\bar 8$

拟随机数的优势在于它的收敛速度更快。如果把估计误差写成计算次数n的函数的话，使用拟随机数会以更高阶的速度收敛。以下是采用拟随机数的误差收敛函数，d表示将要生成的随机数的维度：

1型蒙特卡洛问题（拟随机数）

E r r o r \propto {(\ln n)}^{d} n^{- 1}

2型蒙特卡洛问题（拟随机数）

E r r o r \propto n^{- \frac{1}{2} - \frac{1}{2 d}}

离散和连续布朗运动

随机游走

$W_t$ $W_0=0$ $s_i$ $s_i$ $s_i$ 只能是+1或-1，且概率均为1/2。离散布朗运动定位为按如下公式生成的随机游走序列

W_{t} = W_{0} + s_{1} + s_{2} + c d o t s + s_{t}

其中t=0，1，2，...。

$W_t$ $W_0$ $W_1$ $W_2$ $W_t=s_0+s_1+\cdots+s_t$ 的期望为0，方差为t。

很多随机游走模型的应用都是关于逃逸时间，又称为首次到达时间。设a，b为正整数，一个初始位置为0的随机游走序列，首次达到[-b,a]区间边缘的时刻，就称为逃逸时间。理论证明在a处（而不是-b处）逃逸的概率为b/(a+b)。对于从区间[-b,a]逃逸所需的时间，其期望值为ab。

连续布朗运动

标准随机游走过程在时刻t的期望为0，方差为t。假设在每个单位时间内都可以进行两次位移，且每次位移的时间为1/2个单位时间，则在时刻t的期望仍然会是0，但是方差会变为

V (W_{t}) = V (s_{1} + \dots + s_{2 t}) = 2 t

$1/\sqrt k$ ，这样才能保证方差不变。

$W_t^k$ $\sqrt k$ 的等概率分布，则期望仍然为0，方差为

V ({W_{t}}^{k}) = \sum_{i = 1}^{k t} V (s_{i}^{k}) = \sum_{i = 1}^{k t} [\frac{1}{2} {(\frac{1}{\sqrt{k}})}^{2} + \frac{1}{2} {(- \frac{1}{\sqrt{k}})}^{2}] = k t \frac{1}{k} = t

$W_t^{\infty}$ $B_t=W_t^{\infty}$ 是在t≥0上的随机变量，有三个重要特性

$B_t$ 的均值为0，方差为t；

$t_1<t_2$ $B_{t_2}-B_{t_1}$ $B_{t_1}$ $B_s$ $t_1$ ；

$B_t$ 是一条连续路径。

随机微分方程

有噪声的微分方程

定义 $x_t$ 称为时域连续随机过程。

考虑如下的有噪声的微分方程（SDE）问题：

{\begin{matrix} d y = r d t + σ d B_{t} \\ y (0) = 0 \end{matrix}

$y(t)=rt+\sigma B_t$ 。

SDE是以微分的形式给出的，而在ODE中使用的是导数形式。这是因为像包括布朗运动在内的很多随机过程，它们是连续但不可微的。因此，SDE

d y = f (t, y) d t + g (t, y) d B_{t}

等价于积分形式的表示

y (t) = y (0) + \int_{0}^{t} f (s, y) d s + \int_{0}^{t} g (s, y) d B_{s}

公式中第二个积分项称为Ito积分。

$a=t_0<t_1<\cdots<t_{n-1}<t_n=b$ 是在区间[a,b]上的网格点，黎曼积分定义为如下的极限

\int_{a}^{b} f (t) d t = lim_{Δ t \to 0} \sum_{i = 1}^{n} f ({t_{i}}^{'}) Δ t_{i}

$\Delta {t_i} = {t_i} - {t_{i - 1}},{t_{i - 1}} \le {t_i}' \le {t_i}$ 。类似地，Ito积分定义为如下的极限：

\int_{a}^{b} f (t) d B_{t} = lim_{Δ t \to 0} \sum_{i = 1}^{n} f (t_{i - 1}) Δ B_{i}

$\Delta {B_i} = {B_{{t_i}}} - {B_{{t_{i - 1}}}}$ ，是布朗运动中单次位移的距离。

$B_t$ $dB_t$ 称为白噪声。

Ito公式

若y=f(t,x)，则
$d y = \frac{\partial f}{\partial t} (t, x) d t + \frac{\partial f}{\partial x} (t, x) d x + \frac{1}{2} \frac{\partial^{2} f}{\partial x^{2}} (t, x) d x d x$

有许多数值方法用于求解SDE，有Euier-Maruyania方法，Milstein方法，一阶随机龙格-库塔方法，不过多赘述。

三角插值和FFT

DFT和快速傅里叶变换的内容不过多赘述。

三角插值

DFT插值定理

$\Delta t=(d-c)/n$ $t_j=c+j\Delta t$ $x_j$ $x_j$ $t_j$ 。

$F_n$ x是x的傅里叶变换，由于x是y的逆傅里叶变换，x可以表示为

\begin{aligned} x_{j} = \frac{1}{\sqrt{n}} \sum_{k = 0}^{n - 1} y_{k} {(w^{- k})}^{j} = \frac{1}{\sqrt{n}} \sum_{k = 0}^{n - 1} y_{k} e^{i 2 π k j / n} \\ = \sum_{k = 0}^{n - 1} y_{k} \frac{e^{i 2 π k j / n}}{\sqrt{n}} = \sum_{k = 0}^{n - 1} y_{k} \frac{e^{\frac{i 2 π k (t_{j} - c)}{d - c}}}{\sqrt{n}} \end{aligned}

$(t_j,x_j)$ $y_k$ $y_k=F_nx$ $F_n$ $\{x_j\}$ 变换为插值系数。

$y_k=a_k+ib_k$ ，可以将插值函数写成

Q (t) = \frac{1}{\sqrt{n}} \sum_{k = 0}^{n - 1} (a_{k} + i b_{k}) (\cos \frac{2 π k (t - c)}{d - c} + i \sin \frac{2 π k (t - c)}{d - c})

$x_j$ $x_j$ ，实部为

P (t) = P_{n} (t) = \frac{1}{\sqrt{n}} \sum_{k = 0}^{n - 1} (a_{k} \cos \frac{2 π k (t - c)}{d - c} - b_{k} \sin \frac{2 π k (t - c)}{d - c})

$P_n$ 称为n阶三角函数。

当n为偶数时，有

P_{n} (t) = \frac{a_{0}}{\sqrt{n}} + \frac{2}{\sqrt{n}} \sum_{k = 1}^{n / 2 - 1} (a_{k} \cos \frac{2 π k (t - c)}{d - c} - b_{k} \sin \frac{2 π k (t - c)}{d - c}) + \frac{a_{n / 2}}{\sqrt{n}} \cos \frac{n π (t - c)}{d - c}

$P_n(t_j)=x_j$ ，j=0，...，n-1。

当n为奇数时，有

P_{n} (t) = \frac{a_{0}}{\sqrt{n}} + \frac{2}{\sqrt{n}} \sum_{k = 1}^{(n - 1) / 2} (a_{k} \cos \frac{2 π k (t - c)}{d - c} - b_{k} \sin \frac{2 π k (t - c)}{d - c})

FFT和信号处理

使用三角函数计算最小二乘近似问题

正交性和插值

$F_n^{-1}=\bar{F_n}^T=\bar{F_n}$ $\bar F$ $F_n$ 是酉阵，可以给出一种特殊形式的正交阵，可以用来进行插值计算。

定理 $f_0(t)$ $f_{n-1}(t)$ $t_0$ $t_{n-1}$ 为实数。设n×n矩阵

A = [\begin{matrix} f_{0} (t_{0}) & f_{0} (t_{1}) & \dots & f_{0} (t_{n - 1}) \\ f_{1} (t_{0}) & f_{1} (t_{1}) & \dots & f_{1} (t_{n - 1}) \\ ⋮ & ⋮ & ⋮ \\ f_{n - 1} (t_{0}) & f_{n - 1} (t_{1}) & \dots & f_{n - 1} (t_{n - 1}) \end{matrix}]

为实数正交阵，对y=Ax，以下函数

F (t) = \sum_{k = 0}^{n - 1} y_{k} f_{k} (t)

$(t_0,x_0)$ $(t_{n-1},x_{n-1})$ $F(t_j)=x_j$ ，j=0，...，n-1。

令y=Ax，可以直接给出插值函数

\begin{aligned} F (t) = \frac{1}{\sqrt{n}} y_{0} + \sqrt{\frac{2}{n}} y_{1} \cos \frac{2 π (t - c)}{d - c} + \sqrt{\frac{2}{n}} y_{2} \sin \frac{2 π (t - c)}{d - c} \\ + \sqrt{\frac{2}{n}} y_{3} \cos \frac{4 π (t - c)}{d - c} + \sqrt{\frac{2}{n}} y_{4} \sin \frac{4 π (t - c)}{d - c} + \dots + \sqrt{\frac{1}{n}} y_{n - 1} \cos \frac{n π (t - c)}{d - c} \end{aligned}

用三角函数进行最小二乘拟合

DFT方法可以对区间[0,1]上的均匀采样点进行三角插值，假设n为偶数

P_{n} (t) = \frac{a_{0}}{\sqrt{n}} + \frac{2}{\sqrt{n}} \sum_{k = 1}^{n / 2 - 1} (a_{k} \cos 2 π k t - b_{k} \sin 2 π k t) + \frac{a_{n / 2}}{\sqrt{n}} \cos n π t

公式中的项数为n，等于采样点的数目。采样点越多，需要用到的sin和cos函数项就越多。

$a_k$ $b_k$ 都是线性的，所以可以用法线方程求解。

$x_j$ $x_j$ $t_j=j/n$ $f_0(t)$ $f_{m-1}(t)$ ，拟合的结果如下式

P_{m} (t) = \sum_{k = 0}^{m - 1} c_{k} f_{k} (t)

$P_m$ 在最小平方误差的标准下拟合采样点。

$c_0$ $c_{m-1}$ ，使得方程

\sum_{k = 0}^{m - 1} c_{k} f_{k} (t_{j}) = x_{j}

的误差尽可能小，用矩阵形式表示，为

A_{m}^{T} c = x

$A_m$ $A_m^T$ 的各列是相互正交的，对于变量c的法线方程

A_{m} A_{m}^{T} c = A_{m} x

${A_m}A_m^T$ 是单位阵。因此，最小二乘的解为

c = A_{m} x

n个点的插值多项式可以表示为

F_{n} (t) = \sum_{k = 0}^{n - 1} y_{k} f_{k} (t)

则解可以表示为

F_{m} (t) = \sum_{k = 0}^{m - 1} y_{k} f_{k} (t)

这个结论说明，对给定的n个数据点，用m<n个三角函数进行最小二乘拟合，只需计算n阶插值，再取前m个项即可。即x的插值系数Ax，在丢弃高频分量后，仍然是x的最佳拟合函数。n阶展开里的前m项，就是使用m个低阶频率可以达到的最佳拟合。这个性质反映了基函数的“正交性”。

推论 $x_0$ $x_{n-1}$ $t_j=c+j(d-c)/n$ $a_0$ $a_1$ $b_1$ $a_2$ $b_2$ $a_{n/2-1}$ $b_{n/2-1}$ $a_{n/2}$ $F_nx$ 为x的插值系数，且

P_{n} (t) = \frac{a_{0}}{\sqrt{n}} + \frac{2}{\sqrt{n}} \sum_{k = 1}^{n / 2 - 1} (a_{k} \cos \frac{2 π k (t - c)}{d - c} - b_{k} \sin \frac{2 π k (t - c)}{d - c}) + \frac{a_{n / 2}}{\sqrt{n}} \cos \frac{n π (t - c)}{d - c}

$P_m$ $a_{m/2}$ $P_n$ 中这一项没有2，这与原书423页的引理10.10有关）

P_{m} (t) = \frac{a_{0}}{\sqrt{n}} + \frac{2}{\sqrt{n}} \sum_{k = 1}^{m / 2 - 1} (a_{k} \cos \frac{2 π k (t - c)}{d - c} - b_{k} \sin \frac{2 π k (t - c)}{d - c}) + \frac{2 a_{m / 2}}{\sqrt{n}} \cos \frac{n π (t - c)}{d - c}

$(t_j,x_j)$ ，j=0，...，n-1的m阶最优最小二乘拟合。

声音、噪声和滤波

滤波的用法有两种。第一种是用另一个简单的函数尽可能地拟合原始音频．这是信号的压缩，我们可以只储存m个低频分量。滤波的另一个重要应用是去噪，在一个音乐文件里，音乐和语音中可能会混有高频噪声，需要消除这些高频分量来提高音频的质量。

维纳滤波

$c^Tr$ 的期望应为0。

在一个典型应用中，我们需要从含噪声的信号x中恢复出c。信号c可能是一个重要的系统变量，但是在噪声环境下进行测量。或者如下面例子，c可能是音频采样，我们想从中去除噪声。在20世纪中期Norbert Wiener建议，以最小平方误差为目标，寻找一种最优滤波器来去掉x中的噪声。他建议寻找一个对角阵Φ，使得公式

F^{- 1} Φ F x - c

的模尽可能小，其中F表示DFT变换。它的核心思想是通过傅里叶变换后，在频域进行乘以Φ，再进行逆傅里叶变换。这又称为频域滤波，因为我们处理的是经过傅里叶变换的信号x而不是x本身。为了找到最佳的对角阵Φ，注意到

\begin{aligned} {‖ F^{- 1} Φ F x - c ‖}_{2} = {‖ Φ F x - F c ‖}_{2} = {‖ Φ F (c + r) - F c ‖}_{2} \\ = {‖ (Φ - I) C + Φ R ‖}_{2} \end{aligned}

令C=Fc和R=Fr为傅里叶变换，同时注意到根据噪声的定义

{\bar{C}}^{T} R = {\overset{―}{F c}}^{T} F r = c^{T} {\bar{F}}^{T} F r = c^{T} r = 0

因此模可以写成

\begin{aligned} {‖ (Φ - I) C + Φ R ‖}_{2} = {(\overset{―}{(Φ - I) C + Φ R})}^{T} ((Φ - I) C + Φ R) \\ = ({\bar{C}}^{T} (Φ - I) + {\bar{R}}^{T} Φ) ((Φ - I) C + Φ R) \\ \approx {\bar{C}}^{T} {(Φ - I)}^{2} C + {\bar{R}}^{T} Φ^{2} R \\ = \sum_{i = 1}^{n} {(ϕ_{i} - 1)}^{2} {| C_{i} |}^{2} + {ϕ_{i}}^{2} {| R_{i} |}^{2} \end{aligned}

$\phi_i$ $\phi_i$ 求导，得到

2 (ϕ_{i} - 1) {| C_{i} |}^{2} + 2 ϕ_{i} {| R_{i} |}^{2} = 0

$\phi_i$

ϕ_{i} = \frac{{| C_{i} |}^{2}}{{| C_{i} |}^{2} + {| R_{i} |}^{2}}

但一般情况下我们不知道C或R，所以需要某种近似方法，令X=Fx为傅里叶变换，沿用关于信号与噪声的独立性假设，近似地有

{| X_{i} |}^{2} \approx {| C_{i} |}^{2} + {| R_{i} |}^{2}

可以把最优参数选为

ϕ_{i} \approx \frac{{| X_{i} |}^{2} - {| R_{i} |}^{2}}{{| X_{i} |}^{2}}

$|R_i|^2$ $(p\sigma)^2$ $\sigma$ 为噪声的标准差，p为可调节的接近1的参数。

压缩

离散余弦变换

一维DCT

设n为正整数，一维的n阶DCT变换定义为n×n矩阵C，其元素为

C_{i j} = \frac{\sqrt{2}}{\sqrt{n}} a_{i} \cos \frac{i (2 j + 1) π}{2 n}

其中i，j=0，...，n-1，

a_{i} = {\begin{matrix} 1 / \sqrt{2} i f i = 0 \\ 1 i f i = 1, \dots, n - 1 \end{matrix}

定义 $[x_0,...,x_{n-1}]^T$ $y=[y_0,...,y_{n-1}]^T$ ，其中

y = C x

C = \sqrt{\frac{2}{n}} [\begin{matrix} \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} & \dots & \frac{1}{\sqrt{2}} \\ \cos \frac{π}{2 n} & \cos \frac{3 π}{2 n} & \dots & \cos \frac{(2 n - 1) π}{2 n} \\ \cos \frac{2 π}{2 n} & \cos \frac{6 π}{2 n} & \dots & \cos \frac{2 (2 n - 1) π}{2 n} \\ ⋮ & ⋮ & ⋮ \\ \cos \frac{(n - 1) π}{2 n} & \cos \frac{(n - 1) 3 π}{2 n} & \dots & \cos \frac{(n - 1) (2 n - 1) π}{2 n} \end{matrix}]

注意到C是实数正交矩阵，即转置与逆矩阵相同。

定理（DCT插值定理） $x=[x_0,...,x_{n-1}]^T$ $y=[y_0,...,y_{n-1}]^T=Cx$ ，其中C是n阶DCT变换矩阵，则实数函数

P_{n} (t) = \frac{1}{\sqrt{n}} y_{0} + \sqrt{\frac{2}{n}} \sum_{k = 1}^{n - 1} y_{k} \cos \frac{k (2 t + 1) π}{2 n}

$P_n(j)=x_j$ ，j=0，...，n-1。

举例：用DCT变换计算以下点的插值，(0,1)，(1,0)，(2,-1)，(3,0)。

$[1,0,-1,0]^T$ ，计算Cx得到y即为插值函数的系数，计算得到的y为

y = [\begin{matrix} 0.0000 \\ 0.9239 \\ 1.0000 \\ - 0.3827 \end{matrix}]

所以n=4时得到的插值函数为

P_{4} (t) = \frac{1}{\sqrt{2}} [0.9239 \cos \frac{(2 t + 1) π}{8} + \cos \frac{2 (2 t + 1) π}{8} - 0.3827 \cos \frac{3 (2 t + 1) π}{8}]

DCT变换和最小二乘近似

定理（DCT最小二乘近似定理） $x=[x_0,...,x_{n-1}]^T$ $y=[y_0,...,y_{n-1}]^T=Cx$ $y_0$ $y_{m-1}$ ，得到的函数

P_{m} (t) = \frac{1}{\sqrt{n}} y_{0} + \sqrt{\frac{2}{n}} \sum_{k = 1}^{m - 1} y_{k} \cos \frac{k (2 t + 1) π}{2 n}

$\sum\limits_{j = 0}^{n - 1} {{{\left( {{P_m}(j) - {x_j}} \right)}^2}}$ 最小化。

二维DCT和图像压缩

二维DCT

$n^2$ $(s_i,t_j,x_{ij})$ ，其中i，j=0，...，n-1。从最小二乘的观点来看，二维DCT用最优的方法实现了以上的拟合。

$x_{ij}$ $CX^T$ $t_i$ ，对t方向进行一维DCT变换需要按每行来进行；因此再次转置并乘以C，得到

C (C X^{T})^{T} = C X C^{T}

定义 $Y=CXC^T$ 。

定义 $X=C^TYC$ 。

$C^{-1}=C^T$ $X=C^TYC$ $x_{ij}$ 实际上是余弦函数的乘积。

$X = {C^T}YC$ ，2D-DCT插值如下

\begin{aligned} P_{n} (i, j) = x_{i j} = \sum_{k = 0}^{n - 1} \sum_{t = 0}^{n - 1} C_{i k}^{T} y_{k l} C_{l j} = \sum_{k = 0}^{n - 1} \sum_{t = 0}^{n - 1} C_{k i} y_{k l} C_{l j} \\ = \frac{2}{n} \sum_{k = 0}^{n - 1} \sum_{t = 0}^{n - 1} y_{k l} a_{k} a_{l} \cos \frac{k (2 i + 1) π}{2 n} \cos \frac{l (2 j + 1) π}{2 n} \end{aligned}

$P_n(i,j)=x_{ij}$ ，其中i,j=0，...，n-1。

图像压缩

一般对8×8的图像块进行DCT变换，将图像中的能量集中到图像块的左上角。

量化

模q量化

$z = round\left( {{y \over q}} \right)$
$\bar y=qz$

线性量化定义为矩阵

q_{k l} = 8 p (k + l + 1)

其中0≤k，l≤7，其中常量p称为损失参数。

RGB -> GRAY

X_{g r a y} = 0.2126 R + 0.7152 G + 0.0722 B

RGB -> YUV

\begin{matrix} Y = 0.299 R + 0.587 G + 0.114 B \\ U = B - Y \\ V = R - Y \end{matrix}

霍夫曼编码

不过多赘述

改进的DCT和音频压缩

人的听觉系统对频率非常敏感，在压缩和解压过程中产生的任何瑕疵都能被听出来，基于这个原因，音频压缩中通常会用到一些复杂的技巧来掩盖压缩带来的影响。

改进的DCT

定义 $x=(x_0,...,x_{n-1})^T$ ，其变换为

y = E x

其中E为n×n矩阵

E_{i j} = \sqrt{\frac{2}{n}} \cos \frac{(i + \frac{1}{2}) (j + \frac{1}{2}) π}{n}

引理 $c_j$ 为（扩展的）DCT4矩阵的第j列，则有

$c_j=c_{-1-j}$ （各列关于j=-1/2对称），

$c_j=-c_{2n-1-j}$ （各列关于j=n-1/2反对称）。

$c_{n\over 2}$ $c_{{5\over 2}n-1}$ $c_j$ $c_i$ ，如下图所示。

定义 $x=(x_0,\cdots,x_{2n-1})^T$ 的修正的DCT（MDCT）为n维向量

y = M x

其中M是n×2n矩阵

M_{i j} = \sqrt{\frac{2}{n}} \cos \frac{(i + \frac{1}{2}) (j + \frac{n}{2} + \frac{1}{2}) π}{n}

其中，0≤i≤n-1，0≤j≤2n-1。

与前面讲到的DCT变换的最大区别是： MDCT把2n维的向量转换为了n维的向量。基于这个原因， MDCT不是可逆的，但是通过重叠这些长度为2n的向量，也可以达到可逆的目的。

我们可以把MDCT的变换矩阵M用DCT4的列来表示，然后用上面的引理来简化，得到

\begin{aligned} M = [c_{\frac{n}{2}} \dots c_{\frac{5}{2} n - 1}] = [c_{\frac{n}{2}} \dots c_{n - 1} | c_{n} \dots c_{\frac{3 n}{2} - 1} | c_{\frac{3 n}{2}} \dots c_{2 n - 1} | c_{2 n} \dots c_{\frac{5}{2} n - 1}] \\ = [c_{\frac{n}{2}} \dots c_{n - 1} | - c_{n - 1} \dots - c_{\frac{n}{2}} | - c_{\frac{n}{2} - 1} \dots - c_{0} | - c_{0} \dots - c_{\frac{n}{2} - 1}] \end{aligned}

对n=4的MDCT矩阵

M = [c_{2} c_{3} | c_{4} c_{5} | c_{6} c_{7} | c_{8} c_{9}] = [c_{2} c_{3} | - c_{3} - c_{2} | - c_{1} - c_{0} | - c_{0} - c_{1}]

令A和B表示DCT4矩阵的左半部分和右半部分，则E=[A|B]。定义一个置换矩阵R，可以把矩阵的各列顺序颠倒：

R = [\begin{matrix} 1 \\ . . . \\ 1 \end{matrix}]

当一个矩阵右乘以矩阵R，则它的各列的顺序会颠倒。若矩阵左乘以R，则其各行的顺序会颠倒。

利用置换矩阵，可以进一步简化MDCT矩阵的表示

M = (B | - B R | - A R | - A)

其中AR和BR是列序颠倒后的A和B矩阵。

MDCT变换可以用DCT4来表示，令

x = [\begin{matrix} x_{1} \\ x_{2} \\ x_{3} \\ x_{4} \end{matrix}]

$x_i$ 是一个n/2维的向量（注意n为偶数）

M x = B x_{1} - B x_{2} - A R x_{3} - A x_{4} = [A | B] [\begin{matrix} - R x_{3} - x_{4} \\ x_{1} - R x_{2} \end{matrix}] = E [\begin{matrix} - R x_{3} - x_{4} \\ x_{1} - R x_{2} \end{matrix}]

$Rx_2$ $Rx_3$ $x_2$ $x_3$ 的元素顺序颠倒。这样我们就可以用正交矩阵E来表示M。

由于MDCT变换的矩阵M是一个n×2n矩阵，而不是方阵，因此它是不可逆的。但是，两个相邻的MDCT 的总阶数为2n，因此若把它们组合起来，就可以完美地重建出输入数据x。

$N=M^T$ ，是M的转置。

N_{i j} = \sqrt{\frac{2}{n}} \cos \frac{(j + \frac{1}{2}) (i + \frac{n}{2} + \frac{1}{2}) π}{n}

N = [\begin{matrix} B^{T} \\ - R B^{T} \\ - R A^{T} \\ - A^{T} \end{matrix}]

因为E是正交阵，可以得到

\begin{aligned} A^{T} A = I \\ B^{T} B = I \\ A^{T} B = B^{T} A = 0 \\ R^{2} = I \end{aligned}

N M [\begin{matrix} x_{1} \\ x_{2} \\ x_{3} \\ x_{4} \end{matrix}] = [\begin{matrix} B^{T} \\ - R B^{T} \\ - R A^{T} \\ - A^{T} \end{matrix}] [B x_{1} - B x_{2} - A R x_{3} - A x_{4}] = [\begin{matrix} x_{1} - R x_{2} \\ - R x_{1} + x_{2} \\ x_{3} + R x_{4} \\ R x_{3} + x_{4} \end{matrix}]

在音频压缩中，MDCT变换的向量对应于重叠的数据段。压缩误差会使每段向量之间的连接处产生不连续，而由于向量的长度是固定的，因此这种误差会带来固定频率的噪声。而听觉系统对周期性错误的敏感度要高于视觉系统；实际上，固定周期的错误就是那个频率上的一个声调，而人耳正是感知的各种声调的。如果数据表示采用重叠的方式，令n为正偶数，

Z_{1} = [\begin{matrix} x_{1} \\ x_{2} \\ x_{3} \\ x_{4} \end{matrix}], Z_{2} = [\begin{matrix} x_{3} \\ x_{4} \\ x_{5} \\ x_{6} \end{matrix}]

$x_i$ $Z_1$ $Z_2$ 有一半长度是重叠的，

N M Z_{1} = [\begin{matrix} x_{1} - R x_{2} \\ - R x_{1} + x_{2} \\ x_{3} + R x_{4} \\ R x_{3} + x_{4} \end{matrix}], N M Z_{2} = [\begin{matrix} x_{3} - R x_{4} \\ - R x_{3} + x_{4} \\ x_{5} + R x_{6} \\ R x_{5} + x_{6} \end{matrix}]

$NMZ_1$ $NMZ_2$ $x_3$ $x_4$ ]：

[\begin{matrix} x_{3} \\ x_{4} \end{matrix}] = \frac{1}{2} {(N M Z_{1})}_{n \dots 2 n - 1} + \frac{1}{2} {(N M Z_{2})}_{0 \dots n - 1}

位量化

通过对信号的MDCT进行量化，可以实现音频信号的有损压缩。通过推广图像压缩中使用的量化方法，可以更好地控制信号表示所需的位数。

首先确定一个实数开区间(-L，L)，假设我们要用b位来表示(-L，L)上的一个数，而且可以接受一定程度的误差。我们可以用一个位来表示符号，并用b—1位来对数值进行量化。公式为:

(-L,L)的b位量化

$z=round({y\over q})$ $q={2L\over 2^b-1}$
$\bar y = qz$

特征值与奇异值

幂迭代方法

在求解方程中遇到的威尔金森多项式，难以寻找特征值。本节介绍的方法基于对矢量乘上矩阵的高阶幂，这个向量随着幂的升高会变成特征向量，在后面我们将对这个方法进行精化，但是这就是最复杂方法的主要思想。

幂迭代

幂迭代的动机是与矩阵相乘可以将向量推向主特征向量的方向。

定义：令A是一个m×m矩阵，A的占优特征值是量级比矩阵A所有其他特征值都大的特征值λ，如果这样的特征值存在，与λ相关特征向量被称为占优特征向量。

幂迭代本质上是一个在每步上进行归一化的不动点迭代。和FPI一样，该方法线性收敛，意味着在收敛过程中，在每个迭代步骤中误差以常数因子降低。在本节的随后部分，将讨论幂迭代具有二阶收敛的变种，该方法称为瑞利商迭代。

$x^Tx\lambda=x^TAx$ 的解，或者

λ = \frac{x^{T} A x}{x^{T} x}

这就是瑞利商，给定特征向量近似，瑞利商是特征值的最优近似。

幂迭代

$x_0$
for j = 1，2，3，...
$u_{j-1}=x_{j-1}/||x_{j-1}||_2$
$x_j=Au_{j-1}$
$\lambda_j=u_{j-1}^TAu_{j-1}$
end
$u_j=x_j/||x_j||_2$

幂迭代的收敛

幂迭代的逆

幂迭代局限于求解最大（绝对值最大）的特征值，如果幂迭代用于矩阵的逆矩阵，可以找到最小的特征值。

逆向幂迭代

$x_0$ 以及平移s
for j = 1,2,3,...
$u_{j-1}=x_{j-1}/||x_{j-1}||_2$
$(A-sI)x_j=u_{j-1}$
$\lambda_j=u_{j-1}^Tx_j$
end
$u_j=x_j/||x_j||_2$

$(A-sI)^{-1}$ $(A-sI)^{-1}$ $(A-sI)y_{k+1}=x_k$ $\lambda=b^{-1}+s$ 为矩阵A在s附近的特征值．和λ相关的特征向量可由计算直接给出。

瑞利商迭代

瑞利商可以同逆向幂迭代同时使用，我们知道它会收敛到和s最接近的特征值对应的特征向量，如果这个距离很小则收敛速度很快。如果在这个过程中任何步骤里，近似的特征值已知，可以使用该特征值作为s，以加速收敛。使用瑞利商作为逆向幂迭代中更新的平移会得到瑞利商迭代(RQI)方法。

瑞利商迭代

$x_0$
for j=1,2,3,...
$u_{j-1}=x_{j-1}/||x_{j-1}||$
$\lambda_{j-1}=u_{j-1}^TAu_{j-1}$
$(A-\lambda_{j-1}I)x_j=u_{j-1}$
end
$u_j=x_j/||x_j||_2$

QR算法

$R^m$ 空间（见附录A）的一组单位正交基。这激发了同时对m个向量使用幂迭代，其中我们需要保持每个向量都与其他向量正交。

同时迭代

$AI=\bar Q_1R_1$ 或者

[\begin{matrix} A [\begin{matrix} 1 \\ 0 \\ ⋮ \\ 0 \end{matrix}] | & A [\begin{matrix} 0 \\ 1 \\ ⋮ \\ 0 \end{matrix}] | & \dots | & A [\begin{matrix} 0 \\ 0 \\ ⋮ \\ 1 \end{matrix}] | \end{matrix}] = [{\overset{―}{q_{1}}}^{1} | \dots | {\overset{―}{q_{m}}}^{1}] [\begin{matrix} r_{11}^{1} & r_{12}^{1} & \dots & r_{1 m}^{1} \\ r_{22}^{1} & ⋮ \\ ⋱ \\ r_{m m}^{1} \end{matrix}]

$\bar {q_i}^1$ 是幂迭代过程中新的正交单位向量集，其中i=1，...，m。

归一化同时迭代（NSI）

$\bar{Q_0}=I$
for j=1,2,3,...
$A\bar{Q}_j=\bar{Q}_{j+1}R_{j+1}$
end

$\bar {Q}_j$ $r_{11}^j$ $r_{mm}^k$ 是近似的特征值。


xxxxxxxxxx
function [lam, Q] = nsi(A,k)
[m,n] = size(A);
Q=eye(m,m);
for j = 1:k
    [Q,R] = qr(A*Q)     % QR分解
end
lam = diag(Q'*A*Q);      % 瑞利商

无移动的QR算法

$R_j'Q_j$ $A_j$


xxxxxxxxxx
function [lam, Q] = unshiftedqr(A,k)
[m,n] = size(A);
Q=eye(m,m);
Qbar = Q; R=A;
for j = 1:k
    [Q,R] = qr(R*Q)     % QR分解
    Qbar = Qbar * Q;
end
lam = diag(R*Q);         % 瑞利商

实数舒尔形式和QR算法

定义：如果矩阵T为上三角矩阵，且在主对角线上可以出现2×2的块，该矩阵具有实数舒尔形式。

定理 $A=Q^TTQ$ 。

完整的QR算法迭代通过一系列相似变换，将任意矩阵A移动到它的舒尔分解。我们将使用两个步骤，首先将移动的逆向幂迭代思想植入，并加入收缩技术推出移动QR算法，然后推出改进版本允许处理复数特征值。平移版本可以非常直观地写出。每步中都需要使用平移，完成QR分解，然后再平移回来，使用符号表示：

\begin{aligned} A_{0} - s I = Q_{1} R_{1} \\ A_{1} = R_{1} Q_{1} + s I \end{aligned}

注意到

A_{1} - s I = R_{1} Q_{1} = Q_{1}^{T} (A_{0} - s I) Q_{1} = Q_{1}^{T} A_{0} Q_{1} - s I

$A_1$ $A_0$ $A_k$ $A=A_0$ $A_k$ 的右下角元素，这将导致迭代随着矩阵收敛为实数舒尔形式，而将矩阵的最下面一行除了最右边元素之外全部变为0。在这个元素收敛为特征值后，我们去掉矩阵的最后一行和最后一列进行收缩，然后继续寻找其他的特征值。

原书中有平移QR算法和改进版本的matlab代码

上海森伯格形式

定义 $a_{ij}=0$ ，i>j+1，m×n矩阵A是上海森伯格形式。

定理 $A=QBQ^T$ ，并且B是上海森伯格形式。

可以使用用于构造QR分解的豪斯霍尔德反射构造B，但是在这里有一个主要的差别：现在我们关注在矩阵的左侧和右侧乘上反射子H，这是由于最后想得到一个具有相同的特征值的相似矩阵。正因为如此，我们必须逐步地将0放入矩阵A中。

原书中给出了对应的matlab代码

奇异值分解

$R^m$ $u_1$ $u_m$ $v_1$ $v_n$ $s_1\ge\cdots s_n\ge 0$ ，满足

\begin{aligned} A v_{1} = s_{1} u_{1} \\ A v_{2} = s_{2} u_{2} \\ ⋮ \\ A v_{n} = s_{n} u_{n} \end{aligned}

$v_i$ $u_i$ $s_i$ 是矩阵A的奇异值。

举例：找出矩阵A的奇异值和奇异向量

A = [\begin{matrix} 0 & - \frac{1}{2} \\ 3 & 0 \\ 0 & 0 \end{matrix}]

\begin{aligned} A v_{1} = A [\begin{matrix} 1 \\ 0 \end{matrix}] = 3 [\begin{matrix} 0 \\ 1 \\ 0 \end{matrix}] = s_{1} u_{1} \\ A v_{2} = A [\begin{matrix} 0 \\ 1 \end{matrix}] = \frac{1}{2} [\begin{matrix} - 1 \\ 0 \\ 0 \end{matrix}] = s_{2} u_{2} \end{aligned}

$u_i$ $v_i$ $s_i$ ，m×n矩阵A的奇异值分解（SVD）为

A = U S V^{T}

找出一般的SVD

引理 $A^TA$ 的特征值非负。

$A^TA$ $s_1^2$ $s_n^2$ $v_1$ $v_n$ $u_i$ ，其中1≤i≤m：

$s_i\neq0$ $s_iu_i=Av_i$ $u_i$ ；

$s_i=0$ $u_i$ $u_1$ $u_{i-1}$ 正交。

$u_i$ $s_i=0$ $s_i\neq0$ ，有

\begin{aligned} {(s_{i} u_{i})}^{T} s_{j} u_{j} = u_{i}^{T} s_{i}^{T} s_{j} u_{j} = s_{i}^{T} s_{j} u_{i}^{T} u_{j} \\ {(A v_{i})}^{T} A v_{j} = v_{i}^{T} A^{T} A v_{j} = λ v_{i}^{T} v_{j} = 0 \\ s_{i}^{T} s_{j} u_{i}^{T} u_{j} = λ v_{i}^{T} v_{j} = 0 \end{aligned}

$R^m$ $u_1$ $u_m$ $AA^T$ 的正交特征向量。

定理 $R^m$ $v_1$ $v_n$ $R^m$ $u_1$ $u_m$ $s_1$ $s_n$ $Av_i=s_iu_i$ $V=[v_1|\cdots|v_n]$ $A^TA$ $U=[u_1|\cdots|u_m]$ $A^TA$ 的单位正交向量。

举例：找出2×2矩阵A的奇异值和奇异向量

A = [\begin{matrix} 0 & 1 \\ 0 & - 1 \end{matrix}]

A^{T} A = [\begin{matrix} 0 & 0 \\ 0 & 2 \end{matrix}]

$v_1=[0,1]$ $s_1^2=2$ $v_2=[1,0]$ $s_2^2=0$ $\sqrt 2$ $u_1$ 定义为

\sqrt{2} u_{1} = A v_{1} = [\begin{matrix} 1 \\ - 1 \end{matrix}]

$u_1=[1/\sqrt2, -1/\sqrt 2]^T$ $s_2$ $u_2$ $u_1$ $u_2=[1/\sqrt2,1/\sqrt2]^T$ ，SVD为

[\begin{matrix} 0 & 1 \\ 0 & - 1 \end{matrix}] = [\begin{matrix} \sqrt{2} / 2 & \sqrt{2} / 2 \\ - \sqrt{2} / 2 & \sqrt{2} / 2 \end{matrix}] [\begin{matrix} \sqrt{2} & 0 \\ 0 & 0 \end{matrix}] [\begin{matrix} 0 & 1 \\ 1 & 0 \end{matrix}]

注意SVD分解并不唯一，还有另外一种分解

[\begin{matrix} 0 & 1 \\ 0 & - 1 \end{matrix}] = [\begin{matrix} - \sqrt{2} / 2 & \sqrt{2} / 2 \\ \sqrt{2} / 2 & \sqrt{2} / 2 \end{matrix}] [\begin{matrix} \sqrt{2} & 0 \\ 0 & 0 \end{matrix}] [\begin{matrix} 0 & - 1 \\ 1 & 0 \end{matrix}]

对于对称矩阵，只需计算自身的特征值和特征向量即可。

SVD的应用

SVD的性质

$A=USV^T$ 是奇异值分解，m×n矩阵A的秩是线性无关行的个数（也可以是线性无关列的个数）

性质1 $A=USV^T$ 的秩是S中非0元素的个数。

性质2 $\left| {\det (A)} \right| = {s_1} \cdots {s_n}$

性质3 ${A^{ - 1}} = V{S^{ - 1}}{U^T}$

性质4：m×n矩阵A可以写成秩为1的矩阵的和

A = \sum_{i = 1}^{r} s_{i} u_{i} v_{i}^{T}

$u_i$ $v_i$ 分别是U和V的第i列。该性质是SVD的低秩近似性质，对于A的最优最小二乘近似为保留的前p项，p为矩阵A的秩，p≤r。

降维

$a_1$ $a_n$ $a_1$ $a_n$ ，其中新向量的维数p<m，同时最小化该过程中引入的误差。通常我们从一组均值为0的向量开始，如果向量的均值不为0，我们可以减去均值实现，并在后面的阶段再把均值加回去。

$A=[a_1|\cdots|a_n]$ $A=USV^T$ $e_j$ $Ae_j=a_j$ ，使用性质4的秩近似

A \approx A_{p} = \sum_{i = 1}^{p} s_{i} u_{i} v_{i}^{T}

$a_j$ $u_1$ $u_p$ 张成

a_{j} = A e_{j} \approx A_{p} e_{j}

$e_j$ $<u_1,\cdots,u_p>$ $u_1$ $u_p$ $a_1$ $a_n$ $a_i$ $A_p$ $a_1$ $a_n$ $A_p$ 。

压缩

$u_i$ $v_i$ $s_i$ $s_i$ 。每一项在展开中需要2n+1个数字保存或者传输。

计算SVD

如果A是一个实数对称矩阵，SVD可以消减为本章前面讨论的特征值计算问题。在这种情况下单位特征向量构成正交基。如果我们定义矩阵V将单位特征向量保存在列上，则AV=US表示特征向量方程，其中S是对角线矩阵，保存特征值的绝对值。对于U的操作和V一样，但是如果特征值是负，则要变换列的符号，由于U和V是正交矩阵，

A = U S V^{T}

是矩阵A的奇异值分解。

$A^TA$ ，但当A的条件数很大时，该方法可能不适合使用。

考虑如下矩阵

B = [\begin{matrix} 0 & A^{T} \\ A & 0 \end{matrix}]

注意，B是对称的(m+n)×(m+n)矩阵，因而，它具有实数特征值和作为基的特征向量。令[v,w]表示一个(m+n)向量，对应B的特征向量，则

[\begin{matrix} A^{T} w \\ A v \end{matrix}] = [\begin{matrix} 0 & A^{T} \\ A & 0 \end{matrix}] [\begin{matrix} v \\ w \end{matrix}] = λ [\begin{matrix} v \\ w \end{matrix}]

$A^T$ 得到

A^{T} A v = λ A^{T} w = λ^{2} v

$A^TA$ $\lambda^2$ $A^TA$ $A^TA$ 。

因而，另一种更好的计算奇异值和奇异向量的方法是，将对称矩阵B转化为上海森伯格形式。由于对称，上海森伯格形式的矩阵和三对角线矩阵等价，则该方法与平移QR算法一样可以找出特征值，这是奇异值的平方，特征向量中n个最大的项是奇异向量v、尽管这种方式看起来使得矩阵的大小加倍，但是它避免对条件数的不必要提高。

最优化

最优化指的是找出实数函数的极大或者极小值，该函数称为目标函数。由于定位函数f(z)的极大值与找出函数-f(x)的极小值等价，在推导计算方式时仅考虑最小化问题就足够了。

不使用导数的无约束优化

黄金分割搜索

定义：当区间[a,b]上只有一个极大或者极小值，并且f在其他点上严格升高或降低，连续函数f(x)被称为区间[a,b]上的单峰函数。

定理 $g^k(b-a)/2$ $g = (\sqrt 5 - 1)/2 \approx 0.618$ 。

黄金分割搜索

给出单峰函数f，在区间[a,b]上具有极小值
for i = 1,2,3,...
$g=(\sqrt 5-1)/2$
if f(a+(1-g)(b-a)) < f(a+g(b-a))
b = a+g(b-a)
else
a = a+(1-g)(b-a)
end
end
最终区间[a,b]中包含极小值

持续的抛物线插值

持续抛物线插值

从近似极小值r，s，t开始 $x = {{r + s} \over 2} - {{\left( {f(s) - f(r)} \right)\left( {t - r} \right)\left( {t - s} \right)} \over {2\left[ {\left( {s - r} \right)\left( {f(t) - f(s)} \right) - \left( {f(s) - f(r)} \right)\left( {t - s} \right)} \right]}}$ t = s s = r r = x end

x是对极小值的新的近似。在SPI中，新的x可以替换r，s，t中离当前最远，或者最差的一个点，并重复进行该步骤。对于SPI不能保证收敛，这与黄金分割搜索不同。但是，如果收敛的话，由于该方法更好地使用了函数求值的信息，所以速度会更快。

Nelder-Mead搜索

$R^m$ $x_1$ $x_{n+1}$ $y_1$ $y_2$ $y_{n+1}$ $y_h$ $x_h=x_{n+1}$ $x_h$ $\bar x$ $x_r=2\bar x-x_h$ $y_r=f(x_r)$ $y_r$ $y_1<y_2<y_n$ $x_r$ $x_n$ ，使用函数值对顶点排序，并重复前面的步骤。

具体代码见原书P503页

使用导数的无约束优化

牛顿方法

如果函数连续可微，可以对导数求导，则优化问题可以表示为求解根的问题。

$x^*$ $x^*$ $x^*$ 。对f'(x)=0应用牛顿方法得到如下迭代

x_{k + 1} = x_{k} - \frac{f^{'} (x_{k})}{f^{″} (x_{k})}

当遇到一个多变量函数时，可以使用多变量的牛顿方法，设置导数为0，有

\nabla f = 0

其中

\nabla f = [\frac{\partial f}{\partial x_{1}} (x_{1}, \dots, x_{n}), \dots, \frac{\partial f}{\partial x_{n}} (x_{1}, \dots, x_{n})]

表示梯度。

$F(x)=\nabla f(x)$ $x_{k+1}=x_k+v$ $DF(x_k)v=-F(x_k)$ 的解。梯度的雅可比矩阵DF为

H_{f} = D F = [\begin{matrix} \frac{\partial^{2} f}{\partial x_{1} \partial x_{2}} & \dots & \frac{\partial^{2} f}{\partial x_{1} \partial x_{n}} \\ ⋮ & ⋮ \\ \frac{\partial^{2} f}{\partial x_{n} \partial x_{1}} & \dots & \frac{\partial^{2} f}{\partial x_{n} \partial x_{n}} \end{matrix}]

这是函数f的海森矩阵，因而牛顿步骤为

{\begin{matrix} H_{f} (x_{k}) v = - \nabla f (x_{k}) \\ x_{k + 1} = x_{k} + v \end{matrix}

最速下降

$\nabla f$ $-\nabla f$ 就是最速下降方向。

最速下降

for i = 0,1,2,...
$v = \nabla f({x_i})$
$s=s^*$ $f(x-sv)$
$x_{i+1}=x_i-s^*v$
end

共轭梯度搜索

当A为对称正定矩阵时，求解Ax=w等价于找出抛物面的极小值，如在二维中，线性方程组

[\begin{matrix} a & b \\ b & c \end{matrix}] [\begin{matrix} x_{1} \\ x_{2} \end{matrix}] = [\begin{matrix} e \\ f \end{matrix}]

的解是如下抛物面的极小值

f (x_{1}, x_{2}) = \frac{1}{2} a x_{1}^{2} + b x_{1} x_{2} + \frac{1}{2} c x_{2}^{2} - e x_{1} - f x_{2}

这是由于f的梯度可以表示为线性方程组

[\begin{matrix} a & b \\ b & c \end{matrix}] [\begin{matrix} x_{1} \\ x_{2} \end{matrix}] - [\begin{matrix} e \\ f \end{matrix}] = [\begin{matrix} a x_{1} + b x_{2} - e \\ b x_{1} + c x_{2} - f \end{matrix}] = \nabla f = [\begin{matrix} \frac{\partial f}{\partial x_{1}} \\ \frac{\partial f}{\partial x_{2}} \end{matrix}]

在极小值处的梯度为0，正定意味着抛物面的凹面向上。

$- \nabla f(x)$ ，这是函数f在x点的最速下降方向。假设以及选择了一个搜索方向，使用向量d标记。为了在那个方向找到最小值，需要找出α使得函数h(α)=f(x+αd)最小。将导数设置为0以找出最小：

0 = \nabla f \cdot d = (A (x + α d) - {(e, f)}^{T}) \cdot d = (α A d - r)^{T} d

这意味着

α = \frac{r^{T} d}{d^{T} A d} = \frac{r^{T} r}{d^{T} A d}

共轭梯度搜索

$x_0$ $d_0=r_0=-\nabla f$
for i = 1,2,3,...
$\alpha_i=\alpha$ $f(x_{i-1}+\alpha d_{i-1})$ 最小
$x_i = x_{i-1}+\alpha_i d_{i-1}$
${r_i} = - \nabla f({x_i})$
${\beta _i} = {{r_i^T{r_i}} \over {r_{i - 1}^T{r_{i - 1}}}}$
${d_i} = {r_i} + {\beta _i}{d_{i - 1}}$
end

附录

分块乘法

分块计算A×B，分块的唯一要求是：对A矩阵的各列的分组方法要与B矩阵各行的分组方法一致，如